Critical Behaviour of Random Bond Potts Models: A Transfer Matrix Study 
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We study the two-dimensional Potts model on the square lattice in the presence of quenched 
random-bond impurities. For q > 4 the first-order transitions of the pure model are softened due 
to the impurities, and we determine the resulting universality classes by combining transfer matrix 
data with conformal invariance. The magnetic exponent /3/u varies continuously with q, assuming 
non-Ising values for q > 4, whereas the correlation length exponent v is numerically consistent with 
unity. We present evidence for the correctness of a formerly proposed phase diagram, unifying pure, 
percolative and non-trivial random behaviour. 

PACS numbers: 05.70.Jk, 64.60.Ak, 64.60.Fr 
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I. INTRODUCTION 



The effect of quenched bond randomness on a classi- 
cal statistical mechanics system whose pure version un- 
dergoes a second-order phase transition is well under- 
stood. Namely, the so-called Harris criterion states that 
if the critical exponent a pure governing the divergence of 
the specific heat at the transition point of the pure sys- 
tem is negative, weak bond randomness is irrelevant in 
the renormalisation group (RG) sense and the pure fixed 
point (FP) is stable J|. On the other hand, if aP uro > 
the randomness is relevant and causes a cross-over to crit- 
ical behaviour governed by a new random FP nearby, at 
least if the cross-over exponent a pure is small. 

This should be contrasted with the more dramatical 
effects of randomness in the field conjugate to the lo- 
cal magnetisation. Such randomness can eliminate low- 
dimensional phase transitions altogether, and at least it 
always changes the values of the critical exponents [Q. 
For this reason most early research was concentrated on 
field randomness. In this context a particularly popu- 
lar model is the random field Ising model (RFIM) for 
which a classical argument due to Imry and Ma || from 
a simple comparison of the field fluctuations with the sta- 
bilising effect caused by the formation of a domain wall 
concluded, that the lower critical dimension of the RFIM 
is di = 2. 

The issue of quenched bond randomness imposed on a 
system that undergoes a thermal first-order phase tran- 
sition is less studied. An adaptation of the Imry-Ma ar- 
gument can be established by noting that the bond ran- 
domness couples to the local energy density, which differs 
for the two phases that co-exist at the critical point of 
the pure system, in exactly the same way that the ran- 
dom field couples to the local magnetisation in the RFIM. 
Consequently the existence of a non- vanishing latent heat 
for d < 2 can be ruled out. Early work by Imry and Wor- 



tis H furnished a heuristic argument, reminiscent of that 
of the Harris criterion, that the bond randomness indeed 
softens any such phase transition in d = 2 to a continu- 
ous one. A subsequent phenomcnological RG argument 
by Hui and Berker || confirmed that the lower critical 
dimension for random bond tricriticality and end-point 
criticality is d\ = 2. As the dimensionality increases, tri- 
critical points and critical end points emerge from T = 0. 
Finally, a mathematically rigorous theorem by Aizenman 
and Wehr B stated that quite generally for d < 2 an arbi- 
trarily weak amount of quenched bond randomness leads 
to the elimination of any discontinuity in the density of 
the variable conjugate to the fluctuating parameter. 

The question then emerges whether this softening of 
the phase transition can be verified for specific mod- 
els and, if so, what are the universality classes of these 
novel second-order phase transitions. An investigation 
along these lines has recently been initiated by one of us 
[Q , by considering a system of N two-dimensional Ising 
models coupled by their energy operators which, accord- 
ing to mean-field theory (MFT), is supposed to display 
a second-order phase transition. For N > 2, however, 
the RG flow of the model exhibits a runaway behaviour, 
which is characteristic of a fluctuation driven first-order 
transition 0. In this sense the transition is only weakly 
first-order and hence amenable to perturbative calcula- 
tions. On adding weak bond randomness it was found 
that the RG trajectories curl back towards the pure de- 
coupled Ising FP, and consequently Ising exponents are 
expected, up to possible logarithmic corrections. This 
study was extended by Pujol |j] to the case of N coupled 
q-state random bond Potts models for 2 < q < 4, but 
here the universality class of the impurity softened tran- 
sition was found to depend on the coupling between the 
models. 

A more interesting model for studying the effect of 
quenched bond impurities on a first-order transition is 
the g-state random bond Potts model (RBPM) . For q > 4 
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the phase transition of the pure system is first order with 
a latent heat that is an increasing function of q pj. In 
fact, since the transition is first-order already in MFT, on 
the RG level it is controlled by a zero-temperature dis- 
continuity FP with the eigenvalue of the relevant scaling 
operator being y = d Quenched randomness cou- 
pling to the local energy density thus has the eigenvalue 
d — 2(d —y)=d and is strongly relevant, whence an RG 
treatment appears to be problematic. 

The work undertaken until now has therefore mainly 
been numerical. Extensive Monte Carlo (MC) simula- 
tions have been carried out for q = 8 by Chen, Ferrenberg 
and Landau Jll]] confirming the transition softening sce- 
nario outlined above, and finding critical exponents nu- 
merically consistent with those of the pure Ising model. 
Similar conclusions were reached by Domany and Wise- 
man [fiJI for q — 4 and also for the Ashkin- Teller model. 
It thus appears that in a variety of situations the univer- 
sality class of the bond disordered models is that of the 
Ising model, irrespective of the symmetry underlying the 
original model. 



To explain these findings Kardar et al. |13| have pro- 
posed an interface model for the RBPM which, after sev- 
eral approximations, is amenable to an RG treatment 
that is exact on the hierarchical lattice. In the pure model 
the interface exhibits a branching structure with fractal 
dimension at criticality, but when randomness is present 
the critical interface is asymptotically linear. Assuming 
that the vanishing of the interfacial free energy is gov- 
erned by a zero-temperature FP, the Widom exponent fi 
turns out to be independent of q for all sufficiently large 
g, taking the Ising value /it = 1. 

This is in contrast to the perturbative expansion in 
powers of(o__2) investigated by Ludwig and Cardy [|l4j , 
Ludwig [§j|l6) , and Dotsenko et al. JlJ . Using the RG 
approach for the perturbation series around the confor- 
mal field theories representing the pure models, these 
authors find the critical behaviour of the RBPM to be 
controlled by a new random FP which merges with the 
pure FP as q — > 2. Critical exponents are found to de- 
pend continuously on q, at least for (q — 2) small, and 
in the case of the magnetic exponent xh a calculation 
to three loop order yields a prediction which is supposed 
to be very precise even up to q = 3 (T^|. Unfortunately, 
extending these results beyond q = 4 is impossible, even 
in principle, since this is the limiting case in the range of 
minimal conformal theories around which the perturba- 
tive calculations take place. Another interesting impli- 
cation of this line of research is that the local operators 
exhibit multiscaling Jl6| , meaning that correlation func- 
tions of different moments of such operators decay with 
powers that are, in general, independent. 

It has been suggested by Kardar et al. Jl3| and one of us 
that these contrasting theories describe very different 
FPs. Indeed, it can be argued that the interface model 
pertains to the case of strong non-self-dual randomness, 



whilst the (q — 2)-expansion is relevant for weak self-dual 
randomness. Also, even though it may turn out that the 
critical exponents do not depend on q, the central charge 
c evidently must, since even when the critical behaviour 
is controlled by a decoupled Ising FP there is generally 
not just one Ising model but several. 

To resolve this controversy we have undertaken an ex- 
tensive study of the d — 2 RBPM where finite-size data 
obtained from transfer matrix (TM) calculations were 
combined with the powerful techniques of conformal in- 
variance. We have extended the random cluster model 
TMs of Blote and Nightingale |l8| to the case of bond 
randomness whilst taking into account that in the im- 
pure case such TMs do not commute and hence must be 
discussed in terms of their Lyapunov (rather than the 
eigenvalue) spectra. Because of the lack of self-averaging 
the relation between the Lyapunov spectra and the crit- 
ical exponents is inferred indirectly through a cumulant 
expansion, which has the advantage of illustrating the 
multiscaling properties of the correlation functions (ilf | 
explicitly. The number of Potts states q enters our TMs 
only as a continuous parameter, both facilitating the 
comparison with analytical results within the (q — 2) ex- 
pansion and making the interesting regime q > 4 readily 
accessible. 

Although the cumulant expansion yields very appeal- 
ing results in the case of the magnetic exponent it works 
poorly for the thermal one. For reasons yet not fully 
understood such results, when taken at face value, seem 
to hint at a conformal field theory violating the bound 
v > 2/d |lj||2(|. On the other hand, using phenomeno- 
logical RG techniques |2l| we find results consistent both 
with the bound and with the (q — 2)-expansion p5| , p2|. 

Some of our results have been reported in a Letter p3[ , 
where we also described a mapping between the interfa- 
cial models of the RBPM for large q and the RFIM.. 
This mapping, which is asymptotically exact in the limit 
q — > oo, allowed us to establish a schematic phase dia- 
gram for the RBPM unifying pure, non-trivial random, 
and percolative behaviour. In the present paper the ev- 
idence for this phase diagram will be collected and dis- 
cussed. 

The outline of this paper is as follows. In Sect. || we 
define the model and discuss the principles of extract- 
ing physical information from the Lyapunov spectrum of 
the TMs. The proposed phase diagram is reviewed along 
with the translation of the renormalisation group equa- 
tions from the RFIM to the problem at hand. Then, in 
Sect. [II, the TM formalism of Ref. Jl8| is generalised 
to the random case. This relies on the mapping of the 
RBPM to the random cluster model and on two com- 
plementary representations of the connectivity of a row 
of spin in the latter model. By decomposing the TM 
into sparse single-bond TMs we arrive at a highly effi- 
cient algorithm, the implementation of which is consid- 
ered in detail. The magnetic properties can be accessed 
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by adding a ghost site, but we also descibe an alternative 
route in which the two-point correlator is related to a 
disorder operator under duality. The corresponding im- 
plementation of the TMs has a seam spanning the length 
of the cylinder. In the percolation limit the TMs take 
on a particularly simple form, allowing us to obtain very 
accurate results. 

Sect. [Ey| is dedicated to the presentation of our numer- 
ical results. From the scaling of the free energy we find 
evidence that the first-order phase transition is indeed 
softened due to the randomness. The effective central 
charge is determined both at the random FP and in the 
percolation limit. In the latter case we obtain excellent 
agreement with our analytical result. It is shown how 
a cumulant expansion leads to very accurate values of 
xh- These depend continuously on q, and are in per- 
fect agreement with the (q — 2)-expansion at q = 3. For 
larger q the values stay far away from the Ising value, in 
sharp contrast to the results of Ref. |ll| . For q > 8 the 
expansion begins to break down, and we give an argu- 
ment why this must be so in terms of a model of cou- 
pled replicas. The problems encountered when trying to 
extract xt in a similar fashion then lead us to discuss 
the method of extracting physical observables from the 
Lyapunov spectrum in more physical terms. We then 
consider the constraints put on the multiscaling expo- 
nents by a conformal sum rule. A reliable determination 
of xt is furnished by a variant of the phenomenological 
RG scheme, in which the shape of the self-dual surface is 
explicitly taken into account. The criticism of Ref. JToj ] 
recently raised by Pazmandi et al. ]2~i| ] is shown not to 
apply to the RBPM. We conclude the section with a dis- 
cussion of the higher Lyapunov spectrum and its possible 
relation to the (presently unknown) conformal field the- 
ory underlying the model. 

Finally, Sect, [v] contains a discussion of our findings. 
We seek to explain the discrepancy with Ref. Jll| , and we 
discuss other types of randomness relevant to the ques- 
tion whether a first-order pahse transition is softened due 
to impurities. A list of unsettled questions relevant for 
future research is also given. 



II. THE MODEL AND ITS PHASE DIAGRAM 

A. The random bond Potts model 

The g-state Potts model p5| is defined by the reduced 
Hamiltonian 



H = — ^ K,,<\t -t . 

(«■> 



(1) 



where the spins, defined on the vertices of the square 
lattice, can take the values Oi = 1,2, ... ,q, and the sum- 
mation is over all nearest neighbour bonds in the lattice. 



We shall specialise to the ferromagnetic case, where the 
reduced couplings > measure the strength of the 
aligning tendency of nearest neighbour spins. 

Although the free energy of the pure model (ify = K) 
is not known in closed form for general q, a wide range of 
exact results is nevertheless available [Q. In particular 
it is well-known that the model exhibits a second-order 
phase transition for q < 4 and a first-order one for q > 4 
@. 

However, in this paper we are mainly concerned with 
the random bond Potts model (RBPM) for which much 
less is known. Here the couplings Kij are quenched ran- 
dom variables, typically drawn from the symmetric bi- 
nary distribution 



S(K-K 2 )}, 



(2) 



where the ratio between strong and weak bonds R — 
K2/K1 measures the strength of the randomness. For 
the special choice 
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,K 2 



1) 



(3) 



the model is on average self-dual, as discussed in more 



detail in Sect. [II D below. Assuming that the phase tran- 



sition is unique the model is therefore at its critical point 

Other self-dual distributions of the random bonds than 
that of Eq. (||) have also been investigated in order to 
check our results. In particular, we have found the tri- 
nary distribution introduced in Sect. IV B useful, since it 
gives us a clearer idea about the length scale associated 
with the random impurities. 



B. Lyapunov spectrum of the transfer matrix 

The construction of the transfer matrices (TMs) for the 
RBPM is described in detail in Sect. III. It is well-known 
that in the pure case (R 



1) the operator content of the 
conformal field theory (CFT) underlying the model is re- 
lated to the eigenvalue spectrum {Aj(L)}, i = 0, 1, 2, . . ., 
of the TM for a strip of width L through §1 



fi(L)-f (L) 



2nxi 



+ 



(4) 



where fi(L) = — j-lnAj(L) are the generalised free en- 
ergies per site (in units of k^T) and Xi the scaling di- 
mensions of the corresponding operators. Similarly the 
central charge c, measuring the number of bosonic de- 
grees of freedom of the CFT, is related to the finite-size 
corrections to the customary free energy through p9| 



/o(L) = / (co)- — 



(5) 



In the random case the TMs are no longer constant 
but depend on the particular realisation of the random 
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bonds within each row of strip. Accordingly the concept 
of eigenvalues generalises to that of Lyapunov exponents. 
Starting with some suitable initial vector of unit norm 
|vo), the leading Lyapunov exponent can be found by 
the Furstenberg method |3(J 



H 



RFIM 



-J s i s j - X! h^Si - h}^Sj (7) 



A (L) = lim — In 




(6) 



where Tj is the TM acting between rows j 1 — 1 and j. 
The average free energy per site is given as before by 
fo(L) = — 4;Aq(Z). Higher exponents are found by iter- 
ating a set of n vectors {|vfc)}^~g, where a given \v k ) is 
orthogonalised to the set {Iv;)}^ 1 after each multiplica- 
tion by Tj |jl) . Surprisingly, this method works even for 
a non-hermitian TM, and it is numerically shown to be 
independent of the choice of the initial vectors. 

When some symmetry (e.g., spin reversal or duality) 
is manifest in Tj the orthogonalisation can be circum- 
vented by iterating vectors which belong to definite irre- 
ducible components of that symmetry, but the S q permu- 
tational symmetry inherent in the Potts model has been 
lost through the mapping to the random cluste r mod el 



which forms the backbone of our TMs; see Sect. Ill A 



As to the extraction of physical information from the 
spectra, Eq. (|J) is supposed to retain its validity provided 
that c is replaced by the effective central charge c' , that 
in the standard replica formalism is the derivative of c(n) 
with respect to the number of replicas at n = . The 
question to which extent Eq. (Q) also remains valid is by 
no means trivial and we shall dedicate a fair part of the 
subsequent discussion to it. 



C. Phase diagram 

In the limit q —* oo the behaviour of the pure model is 
readily understood (23) . At the self-dual point the parti- 
tion function is dominated by two contributions, namely 
those corresponding to the q completely ordered states 
and the completely disordered state respectively. All 
other configurations are down by powers of l/y/q and 
have recently been enumerated to 10th order in this small 
parameter [g2) . The dominating states have identical free 
energy but different internal energy densities of K and 
for the ordered and the disordered phase respectively, so 
the transition is, as expected, first order. 

The bond randomness is then included through the 
parametrisation e KiJ — 1 — q2+ w *J : where iv™ = ±w and 
w > measures the strength of the randomness. It can 
now be shown |^3| that as q — > oo the model for an in- 
terface between these two phases of the RBPM is exactly 
the same as that of an interface between the spin-up and 
spin-down phases of the RFIM 



with hf- F — ±ft, RF , provided that one translates quanti- 
ties between the two models using the "dictionary" 
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h <-» -wlnq 
2 

J <-> — In q 



h —t In q. 



(8) 



Here t = (T — T c )/T c is the reduced temperature. 

The infinitesimal RG equations can now be inferred 
from the similar results for the RFIM Q. Near d = 2 
they read 

dw/dl = — (rf/2 — l)w + Aw 3 H (9) 

d{\uq)- 1 /dl = -(]n?) _1 ((d - 1) - Aw 2 + ■ • •) (10) 
dt/dl = t{l + Aw 2 + ■■■), (11) 

where A > is a non-universal constant. The RG flows 
for d > 2 and the proposed phase diagram are shown in 
Fig. (see Ref. @ for details). The shaded region of 
non-vanishing latent heat is bounded by the line Rqi of 
tricritical points. This line, controlled by the fixed point 
R at infinite q, merges with the abscissa as d — > 2. At 
q = oo the intcrfacial mapping is exact so that the flows 
along RPi must extend all the way to w = oo, and since 
w^ 1 is known to be a relevant scaling variable at the per- 
colation limit p4|] , Ref. (23) concluded that Rq2 must be 
separated from the percolative behaviour along P1P2 by 
another line of stable FPs emerging from Pi . It was then 
conjectured that this connects on to the line of random 
FPs found in the (q— gi)-expansion [Q. In this paper we 
shall present the evidence for this conjecture for d = 2, 
when qi = 2 and q2 = 4. 



III. THE TRANSFER MATRICES 

In spite of the large amount of high-precision results 
obtained by combining transfer matrix (TM) techniques 
with finite-size scaling for almost any conceivable type 
of pure statistical mechanics system (see, e.g., Ref. ||| 
for a review) the use of TMs in the study of disordered 
systems seems to have attracted rather little interest as 
compared with the complementary approach of Monte 
Carlo simulations. 

A straightforward way of setting up the TMs for the 
q-state Potts model is to use the traditional spin basis 
where the state of a row of L spins is labelled by the q L 
basis states {<Ti, 02, . . . , ox}, f\ = 1, ■ - ■ , <?■ Whilst this 
approach is highly efficient for q = 2,3 it has two major 
shortcomings in the general case. First, the dimension 
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of the matrices becomes forbiddingly for large q, in par- 
ticular making unaccessible the regime of q > 4 which 
is our main concern. Second, the restriction to integer 
values of q is unnecessary and in fact makes it difficult to 
compare numerical results with analytical calculations in 
the (q - 2)-expansion Jl| 0. 

Both these shortcomings can be remedied by writing 
the TMs in the connectivity basis introduced by Blote 
and Nightingale p8| . In this representation the dimen- 
sion of the TMs is independent of q which enters only 
as a continuous parameter. In fact, the number of basis 
states is asymptotically ~ 4 L (or ~ 5 L upon imposition 
of a magnetic field) with a rather small coefficient of pro- 
portionality, in practice making this basis the preferred 
choice for all but the Ising model (q = 2). 

We have generalised these TMs to include quenched 
bond randomness, and also devised an alternative 
method of accessing the magnetic properties through the 
introduction of a seam along the strip. Furthermore, 
in the percolation limit the TMs are found to simplify 
in a manner that makes calculations for rather large 
strip widths feasible. For convenience these results are 
presented along with a review of the relevant parts of 
Rcf. p8| thus making our description of the Potts model 
TMs self-contained. 



A. Mapping to the random cluster model 

Introducing an imaginary 'ghost site' with fixed spin 
o"o = 1 the partition function for the Potts model can be 
written as 

z = Yi ( n cxp^cw, ) n «p ) . 

(12) 

where Yiuj) ^ s the usual product over pairs of nearest 
neighbour sites and each site i has been connected to the 
ghost site with a similar notation. The reduced mag- 
netic field Hi, here taken to be site dependent, now enters 
at the same footing as the reduced exchange couplings 
Kij. It should be pointed out, however, that a random 
coupling to the ghost site is not a true random field, since 
the latter would try to force different sites into different 
Potts states and not just into the particular state of the 
ghost site with a site-dependent probability. To avoid 
any confusion we shall therefore specialise to the case of 
a homogeneous field Hi = if. 

The site variables can now be traded for bond vari- 
ables through the mapping to the random cluster model 
introduced by Kasteleyn and Fortuin |36| ]. In terms of 
the variables uu = e Kij — 1 and v = e 11 ^ 1 we arrive at 



e f n t)( n 

GCCGoCCa \{ij)&G H J \(iO>6G / 

(13) 

where C denotes the set of all nearest neighbour bonds, 
£q the bonds from each of the N Potts spin to the ghost 
site, and l(G U Gq) is the number of independent loops 
on the combined graph G U Go . 

The usual construction of the transfer matrix T for 
a strip of width L seems to be obstructed by the non- 
local factor 1(GUGq), but this can be taken into account 
by choosing a basis containing information about which 
sites of a given row are interconnected through the part of 
the lattice below that row (including connections via the 
ghost site). This leads us to the concept of connectivity 
states, which we consider next. 

B. The connectivity states 

In order to determine the number of loop closures in- 
duced by appending a new row of L sites along with the 
corresponding L connections to the ghost site to the top 
of G U Go, we need information about how the sites in 
the top row of G U Go were previously interconnected. 
This information is comprised in the connectivity state 
(i\i2 ■ ■ ■ «l), where it = if site t is connected to within 
the combined graph G U Go and, otherwise, i r = i s is a 
(non-unique) positive integer if and only if sites r and s 
are connected within G. 

Whilst this 'index representation' is useful for deter- 
mining whether a newly appended bond does or does not 
close a loop, and thus will allow us to explicitly construct 
the single-bond TMs in the next subsection, a one-to-one 
mapping to the set of consecutive integers {1,2, . . .} is 
clearly needed to define a 'number representation' which 
will enable us to label the entries of the TM and thus 
to perform actual computations. These representations 
and the mapping were supplied by Ref. ]lq | as were the 
determination of the number of connectivity states (d^ 
with and cl without a magnetic field). We shall review 
the necessary details and also give details on the con- 
struction of the inverse of the mapping just mentioned. 

Consider first the case of if = where all ghost bonds 
carry zero weight (v = 0). The connectivity states then 
have all it > and can be recursively ordered by not- 
ing that the index representation is well-nested, i.e., for 
r < s < t < u 

(i r = it) A (i 3 = i u ) =S> i s = i t . (14) 

It follows that if we define the cut function p(i\%2 ■ ■ ■ i£) 
to be the smallest t > 1 such that ii — i t , if such a t 
exists, and L + 1 otherwise, the left (12*3 • • • V— 1) an d 
right (i p i p+ i . . . i£) parts of the index representation are 
both well-nested. A complete ordering of the well-nested 
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sequences is now induced by applying the cut function 
first to the whole sequence, then recursively to its right 
and finally to its left part. 

More precisely, the mapping from the index to the 
number representation is effected by 

( 1 if L < 1 

a(ixt2 ■ ■ ■ ih) = < c L , k -i + [a(i k ■ ■ .i L ) - l]c fc _ 2 (15) 
y +c(*2 ■ ■ • ik-i) otherwise, 

where k = p(i\ . . . i L ) and c n j = J2\= 2 c i-2 c n-i+i with 

(2n)\ 

Cn = Cn,n+X = ~r f — TTT (16) 

n\[n + 1)! 

giving the number of well-nested n-point connectivities 
pj|p7| . Explicit values are shown in Table |. 

To consider the general case of v ^ we remark that 
the subsequence of non-zero indices is still well-nested. A 
complete ordering of an index representation (ii%2 ■ ■ ■ ih) 
with precisely s zero indices is then induced by first or- 
dering according to the value of s, then lexicographi- 
cally ordering the zeros, and finally using the ordering 
of the well-nested subsequence (i Pl i P2 ■ ■ ■ i PL -s) given by 
Eq. (|l5|). The lexicographic ordering is carried out by 

( I if L = I or s = L 
ip(iii 2 ...i L ) = < V , (*2«3---ii) if k (17) 
[ (V) +V(»2»3...i£) if *1 = 0, 

and the mapping to the number representation is finally 

r(iii 2 •••«!,) = c?l, s -i + [ip(iih ■ ■ ■ ih) ~ 1]cl- s 

+ a(i Pl i P2 ...i Pli _ s ), (18) 

where d Uy i = J2\=o (™) c ™-« witn 




giving the number of general n-point connectivities. 
Again, explicit values are presented in Table |j} 

To construct the inverse mapping, i.e., the one tak- 
ing us from the number to the index representation, we 
solve r = T(i\i2 . . Al) for the indices {i\%2 ■ ■ - ih) by 
performing the following steps. First, the number of 
zero indices is found as s = max{s|<i£ iS _i < t}. Sec- 
ond, perform a slightly modified integer division by writ- 
ing r — cLl, s -i — Qcl-s + R, where the remainder R 
is restricted to take its values in the interval [l,cj,_ s ]. 
From Eq. ( |f8| ) we infer that ip = Q + I and a = R. 
Third, the position of the first (leftmost) zero index 
is given by iL-k — 0, where k — max{fc|(^) < -ip}. 
This procedure of finding the zero indices is then iter- 
ated with V -> = ip ~ (*) until ^ s "> = 1, and 
the remaining s — s' zeros are filled in from the right: 
is'+i = ■ ■ ■ = i s -i = i s = 0. 



It remains to deduce the subsequence of non-zero in- 
dices by inverting a = cr(i pi i p2 . ..i Pl ) with / = L — s. 
After initialising i pi — pi we proceed by recursion as fol- 
lows. First, choose k = Tom{k\ci t k-i + Cfc_2Q-fe+i > cr}- 
If k < I we have then found a connection: i pi — i Pk . 
This procedure of finding the connections is now iter- 
ated on the left (i P2 , . . . , i Pk _ 1 ) and the right (i Pk , . . . , i Pl ) 
parts of the remaining sequence. If k > 2 the assign- 
ment i P2 — P2 is performed. By (modified) integer 
division we then write a — Ci^—i = Qck-2 + R with 
R G [1, q_fc + i] , and pass over the left part of the sequence 

with a — > cr}i2 = R and / — > iP"2 = k — 2, and the right 
part with a — > cr\^ ht = Q + 1 and / —> l^ ht = I — k + 1. 
The recursion stops when for any sequence < 2. If 
then /( m ) = 2 and the sequence is (i Pa , i Pa+1 ) we perform 
the assignment i Pa+1 = i Pa if cr^ m ' — 1 and i Pa+1 = p a +i 
if o-M = 2. 

Any way of constructing the index representation 
(«ii2 . . - ih) will of course reflect the above-mentioned ar- 
bitrariness as to the actual values of the non-zero indices, 
but the particular procedure just outlined is easily seen 
to ensure that all indices are < L. This invariant is useful 
since then any given site t can be disconnected from the 
rest by assigning i t = L + I. 

C. The single-bond transfer matrices 

The amount of computer time necessary for building 
up a long strip by repeated application of the transfer 
matrix T can be enormously reduced by decomposing 
the latter as a product of sparse matrices, each corre- 
sponding to the addition of a single bond to C. 

Specifically we write T = T°T h T v , where T v = 
72" • • • T 2 v Ti is connecting each of the L spin sites in the 
uppermost row of the strip to a new spin site situated 
vertically above it, and T h = 1 ■ ■ ■ T^^T^ is finishing 
the new row of C by appending horizontal bonds between 
each of the nearest-neighbour dangling ends created by 
T v . The matrix 7j 1 imposes periodic boundary condi- 
tions by interconnecting the newly added spins at sites 
L and I. Finally T° = T£ ■ ■ ■ T 2 a Tf furnishes the bonds 
of Cq from each of the new spin sites to the ghost site. 
Each of these single-bond TMs is implicitly understood 
to depend on the particular realisation of the bond and, 
in the case of the field randomness pertaining to the 
bond in question. 

Upon addition of one single bond the summation over 
graphs in Eq. ( |l3| ) is augmented by a sum over the two 
possible states of this new degree of freedom, viz. the 
bond added to C (Cq) can be either present or absent in 
G (Go). Correspondingly each column of the TM has at 
most two distinct non-zero entries. 

Consider first adding a vertical bond by action of 7J V , 
/ G {1, . . . , L}. If the bond is 'present' any given con- 



G 



nectivity state (iii2 ■ ■ ■ i£) of the L uppermost spin sites 
will be left unchanged. In case of an 'absent' bond site 
I will be disconnected, and the number representation 
of the new connectivity state can be found by assigning 
il = L + 1 and using Eq. (|l8|). Interpreting the fac- 
tor of q N in Eq. ( |l3| ) as an extra factor of q going with 
each vertical bond we see that the non-zero entries in 
corresponding to a column with a given connectivity 
number are a diagonal contribution of and a possibly 
off-diagonal contribution of q. In particular the vertical 
bonds do not induce any loop closures. 

Similarly the TM of a horizontal bond T^\ +1 has a di- 
agonal entry of 1 for each column, corresponding to the 
bond being absent. The other non-zero entry corresponds 
to a present bond, and its value depends on whether a 
loop is being closed or not. Given the connectivity state 
(ii . . . iiii+i ■ ■ - ih) of some column in the TM this is de- 
termined by comparing i\ and ii+i: if they are equal we 
get an additional diagonal contribution of it;,z+i corre- 
sponding to a loop closure, whereas if they are different 
there is an off-diagonal entry with value uij+i/q. In the 
latter case the connectivity number is found by assigning 
the value min{i/,i; + i} to all indices that were formerly 
equal to either ii or ii + \ and applying Eq. (|l8|). (The 
reason why we copy the minimum index is to ensure the 
proper handling of spins connected to the ghost site.) 

Finally, the TM of a ghost bond has the same form 
as in the case of a horizontal bond if we make the sub- 
stitutions uy+i — > vi and ii+\ — ► 0. 



D. Magnetic properties 

It is well known, at least in the case of a pure sys- 
tem, that physically interesting quantities like the cen- 
tral charge c as well as the thermal (xt) and the magnetic 
{xh) scaling dimensions can be extracted from the trans- 
fer matrix spectrum. Consider for the moment the case 
of vanishing magnetic field, H = 0. Since connections to 
the ghost site are then generated with zero weight (y = 0) 
such connections can only be present in any row if they 
were already there in the preceding row. In particular, 
noting that in the numbering of connectivities induced 
by Eq. (|l8|) the non-ghost connectivities precede the oth- 
ers, we see that the TM assumes the following block form 



T = 







T 12 
q-22 



(20) 



where superscript 2 (1) refers to the (non-)ghost connec- 
tivities. 

The largest and the next-largest eigenvalues of T turn 
out to be the largest eigenvalue of block T 11 and T 22 
respectively, and from the corresponding (reduced) free 
energies per site fo(L) = — -^A" (i = 1,2) for a strip 



of width L the magnetic scaling dimension can be found 
from the CFT formula 6§1 



/ 22 (i) - f?{L) = 



2itxh 
L 2 



(21) 



Physically this relation to xh can be understood by not- 
ing that by acting repeatedly with T 22 on some initial 
(row) state |vo) ^ one measures the decay of clusters 
extending back to row 0. This must have the same spa- 
tial dependence as the spin-spin correlation function and 
hence be related to xh [pV8| ■ Analogously T 11 measures 
the decay of two-point correlations between pairs of spins 
being interconnected within the random cluster model. 
This is nothing but the energy-energy correlation in the 
strip geometry, and accordingly we expect that 



2ttxt 



(22) 



We have checked the results for xh by constructing a 
realisation of the TM in the presence of a seam spanning 
the length of the cylinder. Our algorithm also merits at- 
tention on its own right since it improves the asymptotic 
number of basis states necessary for finding /q 2 (L) from 
dh — Cl ~ 5 L (the dimension of T 22 ) to Lc L — LA L . 
In practice, however, with the strip widths L accessible 
using present-day computers the two algorithms perform 
more or less equally fast (see Table [| for a comparison) . 

The well-known duality relation for the Ising model 
partition function without a magnetic field is easily ex- 
tended to the case of the Potts model on a cylinder. For 
v = the partition function of the random cluster model, 
Eq. (|l3|), can be rewritten as 




C(G) 



(23) 



where C(G) is the number of independent clusters on G. 
We first stipulate the duality between two very special 
graphs. Namely, the full graph G — C with partition 
function Zi n n({Uij}) — q]J 
to the empty graph G* 



(ij) u ij i s taken to be dual 



with Z, 



empty 



({<}) 



where the number of dual sites ./V* is fixed by the Euler 
relation. 

Establishing the duality then amounts to ascertaining 
that all other graphs have the same weight relative to 
this reference state as is the case in the dual model. In 
the terminology introduced above, duality states that a 
graph configuration G on the original lattice C is dual to 
a configuration G* on the dual lattice C* in which every 
bond of strength Uij being 'present' in G corresponds to 
the dual bond of strength u*j being 'absent' from G* . 

In particular, removing one bond from the full graph 
(relative weight: l/iiy) must correspond to adding the 
corresponding dual bond to the empty dual graph (rela- 
tive weight: u*j/q), meaning that the bond strengths and 
their duals must obey the relation 
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(24) 



When removing further bonds from G it may happen that 
a new cluster is separated from the rest of the graph, 
yielding an additional factor of q. But such a cluster 
formation corresponds precisely to a loop closure on the 
dual lattice, also giving an extra factor of ql Since all 
graph configurations G can be constructed by succesive 
removals of bonds from the full reference state we have 
thus proven the fundamental duality relation pq] 



z({u tJ }) = q cz*({ u ;A), 



(25) 



Ui 3 1S some constant. 



where C = q~ N * Jl( 

A similar duality relation can be established for the 
spin-spin correlation function. As usual we define the 
local order parameter as m 



M a (r) = h5 CT ( r ),a - ^ ) 



0=1,. 



(26) 



In the high temperature phase all components of the or- 
der parameter vanish, whilst in the ordered (low temper- 
ature) phase the Z q symmetry is spontaneously broken 
and one of the components, say a = 1, has a positive 
expectation value. A simple calculation now shows that 
the correlation function G aa (ri,r 2 ) = (M a (ri)M a (r 2 )) is 
proportional to the probability that the points r\ and r 2 
belong to the same cluster. 

In a cylindrical geometry the graphs with r\ and r 2 , 
taken to be at opposite ends of the cylinder, connected 
correspond to dual graphs where clusters are forbidden 
to wrap around the cylinder. This is equivalent to com- 
puting the dual partition function with twisted boundary 
conditions a — > (er+1) mod q across a seam running from 
n to r 2 . By permuting the Potts spin states the shape of 
this seam can be deformed at will as long as it connects 
r*i and r 2 . Duality thus maps the correlation function 
onto a disorder operator 



(M a { n )M a {r 2 )) = I Y[ exp(-K*S a 



(27) 



z* 



where Z* = Z*({K *}) is the dual partition function with 
periodic boundary conditions. 

The construction of the TM in the presence of a seam 
is facilitated by the following observation: If no cluster is 
allowed to wrap the cylinder, each graph contributing to 
the partition function can be associated with a function 
s(j) of the row number j, such that s(j) = k € {1, . . . , L} 
means that in row j no horizontal bond connecting sites 
k and k + 1 (mod L) is present. For obvious reasons we 
shall refer to s as the virtual seam. We can then write 
the TM in a basis which is the direct product of the L 
possible values of the virtual seam and the costumary cl 
non-ghost connectivities. The virtual seam is initialised 



by assigning to it a definite value in row 0, viz. s(0) = L 
for all graph configurations of that row. 

The single-bond TM of a vertical bond is diagonal in 
s, but a present horizontal bond not inducing a loop clo- 
sure may alter the value of the virtual seam. Let us 



recall from Sect. Ill C that to find the connectivity state 
(ii . . . ijij+i ■ ■ - il) giving the row label of T^ [+1 that cor- 
responds to the off-diagonal entry with value ui^i/q we 
would join the two distinct clusters formerly labelled by 
either i; or ii+i. But such a merger would ruin the in- 
variant stated above, unless we move the virtual seam 
at the same time. On the other hand, if % = ii+i and 
s(j) = I we must explicitly prevent a cluster from wrap- 
ping the cylinder by leaving out that extra diagonal con- 
tribution which would otherwise by implied by the con- 
dition ii = ij.fi. In this case the virtual seam is not 
moved. 

To conclude this section we remark that in the case 
of a planar geometry any rt-point Potts correlation func- 
tion can be mapped to a generalised surface tension by 
duality ||-|o|. 



E. The percolation limit 



In the random bond Potts model the couplings «y > 
are quenched random variables, and the critical point 
can be accessed by drawing them from the symmetric bi- 
nary distribution P(u) = ^[S(u — tti) +6(u — U2)], where 



U\U 2 — q- For details, see Sect. IV. Bond percolation can 
be studied in the limit U\ — > 0, u 2 — > 00 of infinitely weak 
and strong bonds respectively. In this limit considerable 
simplifications occur in the TM, rendering computations 
with rather large strip widths feasible. 

In the percolation limit all single-bond TM have only 
one non-zero entry per column. Recall from Sect. Ill C 



that in the general case there are two such entries of 
which one is diagonal and the other is 'non-trivial'. In 
the case of the strong vertical bonds and the weak hor- 
izontal bonds only the diagonal entries survive, so that 



the matrices Strong 



U2I and T^Ljj = 1 both become 



trivial. On the other hand, a weak vertical bond cor- 
responds to a TM having one non-trivial entry of q per 
column, whilst a strong horizontal bond is represented by 
a TM that is U2 times a non-trivial matrix with entries 
of l's and 1/g's. 

The factors of 1*2 multiplying both T^ rong and Zjrong 
are innocuous albeit infinite, since of the 2L single-bond 
matrices constituting the entire T there will on average 
be L strong ones, hence L factors of u 2 . On the level of 
the specific free energy this amounts to an infinite addi- 
tive constant 



ti\L) = -\nu 2 + fl\L) 



(28) 



independent of the strip width L. In particular, the cen- 
tral charge c can be extracted from the finite quantity 
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As we shall see in Sect. IV this quantity can be found 
by measuring the asymptotic growth of the norm of 

(rij=i ^j 11 ) I v o)j where |v ) is some largely arbitrary ini- 
tial vector. In the percolation limit the TM turn out to 
be so sparse that after a very few iterations the resulting 
vector has only one non-zero component. Computation- 
ally this means that it is sufficient to store the row index 
of that non-zero component as well as its norm. Both 
time and memory requirements are thus enormously re- 
duced, allowing us to access larger system sizes. 

The disadvantage of this projective quality of the per- 
colation point TM is that neither the thermal nor the 
magnetic scaling dimensions can be found from the Lya- 
punov spectrum. In the case of xh an initial vector in 
the T 22 sector will rapidly decay to zero, thus invalidat- 
ing the procedure for finding /g 2 (L), and the alternative 
of using a seam is obstructed by the fact that disallowing 
the entry in the horizontal bond TM that corresponds to 
a cluster wrapping the cylinder is imcompatiblc with the 
argument of pulling out an overall factor of U2 from the 
TM. 



IV. RESULTS 

A. Softening of the transition 

Before attempting to determine the universality classes 
of the RBPM it is essential to make sure that quenched 
bond randomness indeed renders the phase transitions 
second order. For q > 4 the pure system has a first-order 
transition for which the free energy per site is expected 
to scale like 111 



fo(L) = /o(oo) + aL- d cxp(-L/0 



(29) 



where £ is the bulk correlation length and a is an am- 
plitude depending on q. In Fig. ^ we show plots of the 
function 



A(L) = ln[/o(L)-/ (oo)] + dln£ 
~ const — L/£ 



(30) 



for various values of q and the randomness strength R. 
These plots are rather sensitive to the value of /o(oo), but 
although this is only known exactly for the pure model 
JlO[ it can nevertheless be determined with sufficient ac- 
curacy from the parabolic fits described in Sect. IV B 
below. 

For q = 8 the finite correlation length of the pure sys- 
tem (£ ~ 70) is seen to be rendered effectively infinite 
(£ ~ 10 3 ) upon imposition of the randomness, whilst the 
transition of the Ising model (q — 2) simply stays second 
order. Despite the simplicity of these plots we also find a 
fair agreement with the recently found analytical values 



of £ for the pure systems 0]; near q = 4 these assume 
the simple form 



£ ~ — exp 



(31) 



Another criterion for distinguishing between first and 
second-order phase transitions is the values of thef effec- 
tive) exponents xjj and xt as found from Eq. (|2lj) and 
(p2|) respectively. Generally speaking, for pure systems 
with q > 4 these equations give rise to rather poor fits 
which however have extrapolated values of the effective 
exponents that are in the vicinity of, and slightly be- 
low, zero, whereas when randomness is added the fits are 
much better and yield exponents in the interval ]0,2[. 
In view of the problems justifying such fits in the ran- 
dom case (see below) this evidence for a softening of the 
transition is however not to be taken too seriously. 



B. Central charge at the random FP 

The free energy per site /^(L) for the RBPM on long 
strips of width L is readily found from Eq. (||) applied to 
the T 11 sector of the TM. We have performed extensive 
simulations for various values of q and the randomness 
strength R, though in most cases R = 2 was found to 
describe the random FP adequately. 

Representative samples of our data are shown in Table 
||. For each run a normalised initial vector |vo) was pre- 
pared by choosing its components randomly, and after 
discarding the results of the first 2,000 multiplications 
by T^ 1 in order to eliminate transients, data collection 
was made for each 200 iterations until a strip of a total 
length of m = 10 5 had been built up. For q > 2 a total 
of 100 independent runs were made for 1 < L < 8, and 3 
runs for 9 < L < 12, whilst for the Ising model (q = 2) 
we were able to make 100 runs for 1 < L < 13 by using 
the conventional spin basis. Final results and error bars 
were extracted by computing the mean and the standard 
deviation for the totality of patches of length 200. 

It is not a priori obvious that the Lyapunov exponents 
found from Eq. (||) are independent of the norm used. 
The standard norm in both the spin basis and the con- 
nectivity basis is given by the squareroot of the sum of 
the squared components, and these two are of course not 
identical. To impose the spin basis norm on the connec- 
tivity basis each term in the sum must be weighted by a 
factor q c , where C is the number of clusters in the rele- 
vant connectivity state. We have checked the consistency 
of our results by comparing the first few Lyapunov expo- 
nents obtained from imposing the two different norms on 
the connectivity basis, and we find that not only are the 
results identical but there is even a complete agreement 
of the first three significant digits of the error bars. For 
q = 2 we found that the results using the spin basis and 
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the connectivity basis were consistent, but that the error 
bars obtained using the spin basis were slightly smaller. 

Our results for the free energies of the random-bond 
Ising model agree with, and are more precise than, those 
of de Queiroz 1 @. 

Values of the effective central charge d can be ex- 
tracted from Eq. (|5|) by employing various fitting proce- 
dures. In spite of the relatively slow convergence of both 
two-point fits (L, L+l) and straight-line least-squares fits 
against l/L 2 p2| , iterating such fits yields quite good re- 
sults in the pure model. When randomness is added this 
is no longer so, since rather substantial error bars on the 
first estimates prevent us from efficiently iterating the 
fits. 

A better scheme is to include the leading correction to 
the scaling of Eq. (H), which in the pure case has been 
shown numerically to take the form |f§] 



(32) 



One then performs either three-point fits (L, L + l, L + 2) 
or parabolic least-squares fits against l/L 2 , and because 
of the much faster convergence no iteration is needed [E2I . 
Although a correction proportional to l/L 4 , due to the 
operator TT, must necessarily be present in every system 
that is conformally invariant it can of course not be 
guaranteed to be the dominant one in general. 

In Table III the results of parabolic fits including the 
data points for Lo < L < L ma x have been shown as a 
function of Lq. It is seen that Lq must be chosen large 
enough to justify the omission of higher terms in the se- 
ries ([32]), and small enough to minimise error bars. From 
the special cases of the Ising model and of the percola- 
tion point (see Sect. IV C below) we concluded that the 
choice Lq = 3 is optimal. 



Apart from the results shown in Table III we have also 



performed some runs for q — 1.5, finding, as expected 
from the Harris criterion jjj, no difference between the 
results for the pure and the random model. 

In the intermediate regime 2 < q < 4 our results com- 
pare favourably to those of the (q — 2)-expansion, at least 
up to q — 3. On the other hand, it is evident from Fig. |3| 
that the difference between c for the pure model and c' 
for the random one is of the same order of magnitude as 
our error bars, and only near q = 4, where the expansion 
is expected to break down anyway, are our results able 
to distinguish between the two different behaviours. Ex- 
actly at q = 2 the randomness is marginal and logarith- 
mic corrections to the finite-size scaling forms, Eqs. (0) 



and <^), are expected. Whilst this issue has recently at- 
tracted considerable interest in the case of the critical 
exponents [Hjthe corrections to the central charge are 
much weaker ]45| and accordingly our result is consistent 
with that of the pure Ising model. 

In Fig. ^ we have displayed our results for c' as a func- 
tion of log 10 q for selected values of q € [1.5, 64] . We have 
juxtaposed the results for two strengths of the random- 
ness, namely weak randomness (R — 2, closed circles on 
Fig. |^) and strong randomness (R = 10, open circles). 
For small values of q both randomness strenghts give rise 
to the same c', as witnessed by the overlap of the q = 4 
data points. However, for larger q the R = 2 curve flat- 
tens out and grows slower than logarithmically. Sample 
runs show that the same is true for larger values of R, the 
difference being that the range of g-values for which the 
growth is logarithmic is extended as R is increased. This 
is illustrated by the R = 10 curve's staying above but 
very close to, the percolative result ~ log q (see Sect. IV C 



below) for the whole range of q- values shown on the plot. 
Another way to state this is that for fixed q and vary- 
ing R, the quantity d is an increasing function of R that 
eventually reaches a plateau as R becomes large enough. 
It then appears from Fig. [| that for q < 64 the random- 
ness strength R = 10 is sufficient to reach this plateau. 

These findings are interpreted as follows. According to 
the (q — 2)-expansion the randomness strength R* corre- 
sponding to the random FP is an increasing function of q. 
Assuming this FP to persist as we enter the regime q > 4 
(see Fig. [|) we now claim that the monotonicity of R*(q) 
also holds true when the (q — 2)-expansion breaks down. 
From the RG flows given in Fig. [j] we see that any initial 
value of R G]l, 00 [ will eventually flow to the random FP 
as the system is viewed on larger and larger length scales. 
However, if we start out very far from R* the onset of 
the asymptotic scaling given by Eq. (||) may be deferred 
to much larger length scales than the strip widths L nu- 
merically accessible for our TMs. We therefore expect 
poor scaling for strip widths L < L max - Conversely, if 
we choose the strength of the randomness as R ~ R* 
the resulting value of c' is expected to be more or less 
independent of the precise choice of R and equal to the 
true value of the central charge. But in our simulations 
we find that this is precisely accomplished by choosing R 
as an increasing function of q. Further justification for 
this interpretation is found from the phenomenological 
RG treatment in Sect. 



IV F 



below. 



A heuristic argument explaining that the "effective" 
c'(R,q) obtained for small values of R is less than the 



Actually they differ by a constant, since Ref. @ defines the 
Hamiltonian as — Xy(ij) Kij(Ti<jj as opposed to our Eq. (|l|). 
Since SiSj = 2Sa i( y j — 1 there is a free energy difference of 
2~K~, which for R = 2 equals 0.91407. 
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"correct" value c'(R*,q) associated with the random FP 
is readily furnished, at least for large q. Namely, from 
Zamolodchikov's c-theorem we know that there exists 
a function c({K}) of the couplings that decreases along 
the RG flow and equals the central charge at the fixed 
points. As a corollary the curves of constant c are or- 
thogonal to the RG flow. In particular, for large q where 
the RG flow is known from the mapping to the RFIM 
(see Eqs. ^ and ([To|)), it is evident from Fig. |l| that 
c'(R 7 q) is equal to c'(R*,q*) for some q* < q. Since our 
numerical results indicate that c'{R*,q) is an increasing 
function of q the proposition follows. 

To check our results for d we have also made 100 in- 
dependent runs for each of the strip widths 1 < L < 8 
where the random bonds were drawn from the trinary 
distribution 



P(K) 



(l-2p)5(K- 



-S(K 
K*), 



(33) 



where K\ and Ki = lOOOifi satisfy the criterion (Q) and 
(exp K* — l) 2 = q. Here p <C 1 is the strength of the ran- 
domness. Of course this realisation of the randomness 
also preserves self-duality, and hence the model is again 
at its critical point p7| . 

Numerical results for c' using trinary randomness are 



shown in Table [V and they are consistent with the bi- 
nary results given above, again provided that p is in- 
creased as we go to larger and larger q. In particular 
it is reassuring to verify that we seem to probe the true 
random behaviour when 2/p (the length scale associated 
with this randomness) is comparable to the correlation 
length of the pure system (pi]). 

An interesting question is whether the asymptotic 
value of c' is approached from above or below when the 
system is viewed on larger and larger length scales. For 
models exhibiting reflection positivity Zamolodchikov's 
c-theorem |fi"6]j ensures that the convergence is from 
above. In particular the condition of positivity holds true 
for unitary models, whereas for a random model it may 
well fail to be fulfilled. Indeed, in the case of the RBPM 
a perturbative calculation |TJ] suggests that the conver- 
gence may be from below in some cases. 

In order to discuss this point the parabolic fits versus 
1/L 2 employed above are no longer appropriate. Apart 
from speeding up the rate of convergence to a point where 
information about its direction becomes obliterated due 
to error bars the inclusion of higher-order corrections to 
the finite-size scaling form (|J) may have the effect of re- 
versing this direction. E.g., in the case of the pure Ising 
model it is found |l2| that the estimators obtained from 
parabolic fits converge from below, whereas the corre- 
sponding linear fits (i.e., without the 1/L correction) 
yield estimators that converge from above in accordance 
with the theoretical prediction. 

In Table [v] we show the results of such linear least- 



squares fits for several values of q. The randomness 
strength R was chosen in accordance with the considera- 
tions given above. It appears that in all cases the finite- 
size estimators converge towards the asymptotic values 
of Table III from above. 

We remark that values of c' similar to ours have re- 
cently been reported by Picco |^7j . For q = 8 this author 
found d — 1.45 ± 0.06 which agrees with our result of 
respectively c' — 1.52 ± 0.02 for binary randomness of 
strength R = 10, and d = 1.51 ± 0.04 for trinary ran- 
domness of strength p — 0.10. Our observation that 
c' appears to be an increasing function of R, eventu- 
ally reaching a plateau as R becomes large enough, was 
confirmed by Ref. |47j| that used binary randomness of 
strength R = 1Q throughout. Strong evidence was also 
given that c'(q) grows roughly logarithmically with q in 
the regime q € [5,256], but a further discussion of what 
this implies will be deferred to Sect. |v| below. 

It is worthwhile to compare the TM algorithm used in 
Ref. [[l?]] to ours. It was found that the number of distinct 
entries in the pure model TM in the spin basis is 



fe = £££-£*. 



(34) 



82=1 £3 = 1 U = l 



where nij = max(i2, «3, . . . , nii_i), and L designates the 
strip width as usually. Further taking into account the 
2 L different realisations of the binary randomness in each 
strip, recursion relations between the different elements 
of the TM were found by computing a total of (&l) 2 2 l 
polynomials. Since this number of polynomials increases 
rapidly with L high-precision computations could only 
be performed up to L max = 6. The number of iterations 
used to determine fo(L = 6) was similar to ours, whereas 
more iterations were used for the smaller strip widths. 

Evidently this algorithm also has the advantage that 
q enters only as a parameter, thus making accessible any 
value of q for the simulations. However, for large L it 
performs inefficiently, as we will now show. The num- 
bers bh of Eq. ( |34"| ) are by no means unfamiliar. Indeed, 
they are nothing but the total number of L-point connec- 
tivities, including the non- well- nested ones |}{|. Alterna- 
tively they can be viewed as the number of ways that 
L objects can be partitioned into indistinguishable parts 
p0| . With m v parts of v objects each (v — 1,2,.. .) this 
can be rewritten as 



£ 



n 



L\ 



=0 u=l V > 



(35) 



where the primed summation is constrained by the con- 
dition X^^Li vm v = L. From this representation the gen- 
crating function can be immediately inferred 



exp(e* - 1) 



7 , 

0,,,i 

n=0 



n 



(36) 
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Explicit values, found by Taylor expansion of the left- 
hand side, are shown in Table J. Asymptotically the 6l 
are seen to grow faster than L whereas the well-nested 
connectivities only grow as ~ 4 L . 



C. The percolation limit 

In the case of the binary randomness (^|) the perco- 
lation limit is reached by letting (e Kl — 1) — » and 
(c K2 — 1) — > oo whilst maintaining the self-duality cri- 
terion (|3|). The partition function of the random cluster 
model is then dominated by one graph only, viz. the one 
that covers all of the strong bonds and none of the weak 
ones. (Note in particular that the limits R — ► oo and 
q — > oo do not commute.) Expressed in terms of the free 
energy per site this reads 



ppcrc 
JO 



1) 



~N 



(37) 



where B is the number of strong bonds and C is the 
number of clusters in the dominant graph. 

The quenched average over the randomness must be 
taken on the level of the free energy. Evidently, with the 
chosen distribution of the randomness, B = N whence 
the first term is simply a trivial, albeit infinite, constant. 
(Incidentally this is the same constant that was pulled 
out in Eq. (p8|).) On the other hand, the average number 
of percolation clusters is related to a derivative in the 
pure Q-state Potts model || 



C 



d_ 

dQ 



mZ(Q) 



(38) 



2=1 



thus determining the effective central charge c'{q) at per- 
colation as 



c'{q) = lng 



dc(Q) 



dQ 



(39) 



An alternative argument for this relation is furnished by 
the observation that the replicated model is simply the 
Potts model with q n states; differentiating c(q n ) with re- 
spect to the number of replicas n and taking the limit 
n — ► one recovers the result (^9|) . The central charge of 
the pure model is given by an expression due to Kadanoff 



(2-3y)(l + y) 

(2-y) : 



(40) 



where y[Q = 2cos(7ry/2) and < y < 1, and taking the 
appropriate derivative of this we finally arrive at 



c'(q) 



5^/3 

47T 



In q. 



(41) 



As described in Sect. HIE the single-bond TMs in the 
percolation limit have only one non-zero entry per col- 
umn, equal to either q, 1 or l/q. Taken together with 
their projective quality and Eq. (||) for the largest Lya- 
punov exponent it is clear that the free energy, and hence 
the central charge, must be explicitly proportional to In q. 
So it suffices to do the numerics for one value of q ^ 1. 
Because of the simple form of these TMs we were able 
to average /^(L) of Eq. ( p8| ) over 100 strips of length 
m = 10 5 for the range 1 < L < 19. Consequently the 
factor of proportionality could be determined quite ac- 
curately as 0.688 ± 0.003, in excellent agreement with 



5V3 



0.689. 



It is evident from the mapping between bond perco- 
lation and the pure Q = 1 Potts model that the critical 
exponents of the two models are identical: Xt = f and 
Xr = Is" @- Since all correlation functions at perco- 
lation can only take the values and 1, it is also clear 
that different moments of a given correlation function all 
have the same scaling dimension. Thus, in the notation 
of Ludwig [fl6|| , x n = x\ for all n > 1. The pure model 
represents the other trivial extreme case of multiscaling 
behaviour: 



D. The cumulant expansion 

The concept of multiscaling is best understood in terms 
of a simple example p9[ : The random-bond Ising chain. 
From the reduced Hamiltonian Ti = — Y2i=i ^i s i s i+i 
the partition function is easily found as Z = OiLi c ij 
where c-i — 2 cosh Ki . In particular the quenched aver- 
age Z = exp[iV log cl] does not coincide with the most 
exp[Alogc,]. On the other hand, 



probable value Z m . p 
the reduced free energy is F 
F 



DiXi lo g c » so tnat 

m.p. — — -/VlogCj. The free energy is thus self- 
averaging, i.e., it takes on its sample averaged value 
with probability unity in the thermodynamic limit. In- 
deed, by the central limit theorem, F is normally dis- 
tributed, it being a sum of random numbers, whereas 
Z is a product of random numbers and therefore log- 
normally distributed. Similarly the correlation function 
(sisii) — Hili 1 ta nn Ki is non-self-averaging. In partic- 

2 

ular (siSr) 2 ^ (siSr) . 



F n 



In Sect. [II D we related the spin-spin correlation func- 



tion G(m) on a strip of the RBPM to the free energy in 
the presence of a seam of frustrated bonds (or with a 
ghost site). Taking the logarithm of Eq. ( p7| ) and ex- 
ploiting the self-duality of the lattice we have 



Af(L) = / 22 (L) - fl\L) = — hxG(m) 



(42) 



and in the pure system, according to conformal symmetry 
|||], this decays along the strip as 2ttxh/L 2 , cfr. Eq. pi]). 
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When randomness is present A/(L) is a fluctuating quan- 
tity, and since free energies are supposed to be nor- 
mally distributed these fluctuations are 0(l/y/m). Con- 
sequently In G is a self-averaging quantity and G is not 
p9[ , exactly as in the simple example given above. 

In the multiscaling scenario of Ludwig | fl6| different 
moments G(m) n scale with dimensions x n which, as op- 
posed to what is the case in the pure model, are not 
necessarily linear in n. (In this notation xh = x\.) For 
n\ > ri2 we have x ni > x n2 and x ni /ni < x n2 /ri2 (con- 
vexity); pure and percolative behaviour are thus realisa- 
tions of the two possible extremes of multiscaling. 

Since translational invariance is one of the basic as- 
sumptions of conformal symmet ry |43|, the latter only 
refers to the averaged quantities G(m)™ and not to the 
G{m) n themselves. These averages cannot be computed 
directly in a numerical experiment because of the lack 
of self-averaging; this can however be circumvented by 
performing a cumulant expansion 



1 



In G n = nln G + ^n 2 (In G - InG) 2 



(43) 



where each term on the right-hand side is self-averaging 
and can be directly extracted from the statistical fluctu- 
ations in A/(L) between the patches of length 200 into 
which we have divided our stip. 

Quite generally for a stochastic variable X we have 



(expX) = exp^—^kj, 



(44) 



where explicit expressions for the six first cumulants ki 
in terms of the moments m t of X are given by [ |50| 

fci = mi 

k 2 = m 2 — m\ 

&3 = TO3 — 3m2TOi + 2m\ 

ki = rri4 — 4m3mi — 3m?, + \2m2m\ — Qm\ 

&5 = — hm^m\ ~ 1077737772 + 207713m 2 

+ SOmJmi - 60777 2 mJ + 24m' 1 
k e = itiq — 67775m! — 1577747712 + 30m4m 2 — 10m 2 

+ 120m 3 m 2 TOi - 120m 3 mJ + 30m| - 27Qm 2 2 m 2 

+ 360m 2 mi - 120m? 

We have computed these six cumulants of Af(L) for var- 
ious values of R and q, based on 100 independent strips 
of length m = 10 5 and width 1 < L < 7. Sample results 
for R = 2 and q = 3, 8 are shown in Table VI. 

For q = 3 the cumulant expansion converges well. 
The magnitude of the higher cumulants decreases very 
rapidly, especially for L > 3, and reliable estimates for 
the left-hand side of Eq. ( p3| ) can be obtained simply by 
summing the first 3 or 4 cumulants, at least when 77 is 



not too large. Performing parabolic least-squares fits us- 
ing Eq. ( [jl| ) with an l/L 4 correction we thus expect to 
extract quite accurate values of x n at the random FP. 

As q increases the convergence is slower. This is wit- 
nessed by the q — 8 results of Table VI decreasing notice- 



ably slower, both for a definite cumulant as a function of 
L (vertically) and for a definite L as a function of the 
cumulant number (horizontally). The approximation of 
leaving out the higher cumulants in the sum J43| ) thus 
becomes increasingly difficult to justify, and eventually 
the cumulant expansion breaks down. This problem is 
enhanced by the fact that for q > 8 we expect a ran- 
domness strength of R — 2 to be insufficient in order 
to access the true behaviour at the random FP. We are 
thus forced to increase R, whence the fluctuations be- 
come even more violent and the cumulant expansion ac- 
cordingly ill-behaved. 

Our results for x± are shown in Fig. |^. Since error bars 
on the individual cumulants are related to the magnitude 
of the higher cumulants the question of how to assign a 
final error bar to x\ becomes a delicate one. We have 
addressed this issue by averaging the estimates for X\ 
obtained from various parabolic least-squares fits. More 
precisely, the average is calculated from 4 values, namely 
fits with Lq = 3 or 4 and including the first 3 or 4 cumu- 
lants on the right-hand side of Eq. ( fl3| ) . The consistency 
of these 4 values is regarded as a check of the validity of 
the expansion. 

In particular, for q = 3 we find xi(3) = 0.13467 ± 
0.00013 which is 10 standard deviations above the value 



x P ure (3) = § ~ 0.13333 of the pure three-state Potts 
model p3| and at the same time in perfect agreement 
with the result x x (3) = 0.13465 + C(e 4 ) of the (q - 2)- 
expansion @. The Monte Carlo result x x (3) = 0.1337± 
0.0007 of Picco J51J was not able to distinguish convinc- 
ingly between pure and random behaviour. 

For q = 4 our result is xx{4) = 0.1396 ± 0.0005, in nice 
agreement with Picco's preliminary result xi(4) ~ 0.139 
fl52| and decidedly different from the corresponding pure 
value of < lrc (4) = ±. 

As discussed at length in the Introduction a major mo- 
tivation for this work was to determine whether the im- 
purity softened transitions for q > 4 do or do not have 
the critical exponents of the pure Ising model. The data 
of Fig. H clearly show a smooth continuation of the per- 
turbative results |]l6| , |l7| exhibiting no singularity what- 
soever at q = 4 Our result 11.(8) = 0.1415 ± 0.0036 is 
comfortably away from the pure Ising value and provides 
a striking piece of evidence for both our phase diagram 
and the FP structure of the (q — 2)-cxpansion. 

All the results quoted for x\ were computed using 
R = 2. We have checked that other values of R yield 
results consistent herewith, provided that R is chosen 
neither too small, in which case the cross-over length 
ix ~ ex\){\/2Aw 2 ) found from Eq. (^|) becomes too large 
for the random FP to be reached, nor too large, in which 
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case the cumulant expansion breaks down. The same 
holds true when the random bonds are drawn from the 
trinary distribution ( |33| ) with various values for the dilu- 
tion parameter p. 

Because of the positive sign of the second cumulant the 
values of x\ are invariably smaller than those one would 
have obtained without the cumulant expansion (i.e., us- 
ing only the first cumulant). The latter, however, deter- 
mine a universal exponent ckq that describes the asymp- 
totic decay of the spin-spin correlation function in a fixed 
sample at criticality. In terms of the multiscaling expo- 
nents this reads 



a 



dn 



Near q = 2 Ludwig obtained the expansion |T(| 



Id 



(45) 



(46) 



where y is the RG eigenvalue of the energy operator cou- 
pling to the bond randomness. Our results for small frac- 
tional values of (q — 2) are in good agreement with this 
expression, and if one takes into account the logarithmic 
corrections expected exactly at q = 2 it seems that the 
theoretical prediction and the numerical results have a 
common tangent at q = 2. For the physically interesting 
case of q = 3 the agreement is not so good. We believe 
that this apparent discrepancy would be resolved if the 
expansion ([l6]) could be carried through to three-loop 
order as in the case of x\. 

Similar remarks can be made about the higher mo- 
ments of the spin-spin correlation function for which we 
are unable to verify Ludwig's expansion p6[ 



pure 



1 , 
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l)y + 0(y 2 ). 



(47) 



Nevertheless it should be remarked that the fact that the 
higher cumulants do not vanish in itself implies multiscal- 
ing. 

Before concluding this section we should like to give a 
heuristic argument that the cumulant expansion breaks 
down for large q. In a replica formulation we can imag- 
ine the central charge c(n) as a function of the number of 
replicas n. In this notation the central charge of the pure 
and the random systems are c(l) and c'(0) respectively, 
where the prime denotes a differentiation with respect to 
n. The partition function of the replicated strip is then 



Z n = J exp(-nro£/)P(/) df 
Trmc(n) 



= exp ~mLf 



6L 



(48) 



where P(f) is the probability distribution of the free en- 
ergy. Differentiating this expression twice with respect 



to n and taking the replica limit n — > we infer that 
the second cumulant of / contains a term that is propor- 
tional to c"(0). The cumulant expansion is thus expected 
to break down if c(n) has a large curvature at n = 0. 

For 2 < q < 4 the replicas are weakly coupled, since 
c(l) ~ c'(0) @. Hence c"(0) < 1. But when q = 4 + e 
the transition of the pure system goes first order so that 
the function c(n) starts out with slope c'(0) = 1 and 
somehow curves down to assume the value c(l) = 0. Con- 
sequently c"(0) = 0(1) and the higher cumulants begin 
to contribute significantly to the sum (fl3|). Finally, for 
g > 4 we are in the strong coupling regime. We still 
have c(l) = and as our numerical data indicate that 
c'(0) ~ In q it follows that c"(0) > 1. This means that 
the cumulant expansion must break down. 

One may speculate whether the transition actually be- 
comes first-order whenever q n > 4. Clearly this is the 



case for the pure Potts model (l0|, but a similar state- 
ment is true when N Ising models are coupled by their 
local energy density. Namely, in this case an RG anal- 
ysis implies a fluctuation-driven first-order transition 
whenever N > 2, that is to say for 2 N > 4. If this 
conjecture is correct one would then suppose the func- 
tion c(n) to vanish for n > no, where q n ° — 4. Evidently 
such a scenario is in accordance with our observation that 
c"(0) > 1 for q > 4. 



E. The thermal exponent 

Because of the rather striking success of the cumulant 
expansion for x\ one would now expect the thermal ex- 
ponent xt to be similarly related to the fluctuations of 
A./t(L) = — /q-'-(L). Surprisingly, this seems not 

to be the case. Computing the equivalent of oeo, i.e., us- 
ing only the first cumulant, we find the following results 
for different values of q: 0^(2) = 1.028 ±0.001, c#(3) = 
0.91 ±0.01, 0^(4) = 0.81 ±0.02 and c#(8) = 0.65 ±0.01. 
As remarked above the results using more cumulants can 
only be lower. 

This is bad news since the quenched correlation length 
exponent v can be shown quite rigorously to satisfy the 
bound [ppfll 



v > 



d' 



(49) 



or, in our notation, xt > 1- Though the proof of Ref. ]18| ] 
refers to the divergence of the correlation length as the 
critical point is approached, and hence strictly speaking 
does not apply to the system under consideration since 
we work exactly at the critical point, the RBPM is among 
the simplest physical systems for which Eq. (B9h is be- 



lieved to be valid 20 1 . The point is strengthened by not- 
ing that the (q — 2)-expansion yields xt = 1.02 + 0(e 3 ) 
at q = 3 pq]. It is therefore difficult to have confidence 
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in the cumulant expansion for the thermal exponent, and 
independent methods of assessing xt must be devised. 

At this point we note that although the RG equation 
( pd| ) seems to warrant an effective exponent of = 
1 — Aw 2 for q large, this argument is only superficially 
true. Indeed, near q = oo the RG flows must extend 
to infinite w before reaching the random FP, and conse- 
quently an expansion valid for weak randomness is not 
to be trusted. 

The alternative method for finding xt that comes clos- 
est to the spirit of Refs. (l^^^l is that of finite-size scaling 
off the critical point. This is discussed at length in the 
next subsection, and for the moment we concentrate on 
less "obvious" possibilities. 

One of the key points in the construction of the cu- 
mulant expansion was the realisation that the spin-spin 
correlation function was mapped onto a surface tension 
under duality, and hence could be expressed in terms 
of the largest Lyapunov exponent of a TM with twisted 
boundary conditions. Reinterpreting the latter as a free 
energy the self- averaging property was evident, and the 
cumulant expansion correspondingly behaved quite well 
if the fluctuations were not too large. It has recently 
been shown [ |39| , ff0[ that under duality four-point corre- 
lation functions are similarly mapped onto (generalised) 
surface tensions. Presently these duality relations have 
only been worked out for planar graphs, but there is some 
hope that they may be extended to the case of cylindrical 
boundary conditions as well. Taking two of the points as 
nearest neighbours on either end of the cylinder we would 
then recover the energy-energy correlator, and if the cor- 
responding boundary conditions can be implemented in 
the TM xt follows from a cumulant expansion. 

Next, we discuss the method of iterating orthogonal 
vectors in order to extract the second Lyapunov expo- 
nent []3l| in more physical terms. The energy-energy cor- 
relator (Green's function) can be written as 



(b„|b„) (b„|b„) (a„|b„) (b„|a n ) 



(E( ri )E(r 2 )) 



Tr£(n)£(r 2 )exp(-H) 
Tr exp(-W) 



(50) 



Now imagine building up the strip by repeated action 
by the random TMs on some initial state situated at 
r = — oo. When we reach n the system is in a state 
|ao) on which we act with the energy operator to define 
|b ) = E(ri)\&o). After n further iterations these states 
become 



I an) =T n - ■ - T 2 Ti\& ) 
\h n )=T n ---T 2 T l \h Q ). 



(51) 



Defining a new state |b„) by orthogonalising |b n ) with 
respect to |a n ) 



|b„) = \K 



(a n |b r , 



( a n a„ 



(52) 



(a n \8, n 



\ a n , &n 



{E{n)E{r 2 )) - {E( n )) {E(r 2 )). (53) 



Thus the process of orthogonalisation corresponds pre- 
cisely to subtracting off the disconnected part of the cor- 
relation function. 

When n ^> 1 the states |b n ) and | almost iden- 

tical due to contamination and have a huge norm ~ Aq . 
The idea of orthogonalising them is therefore numeri- 
cally extremely unsound. Fortunately a simple calcula- 
tion shows that orthogonalising after n\ iterations and 
then again after n — n\ further iterations is equivalent 
to orthogonalising only once, as above. Hence, by induc- 
tion, we are allowed to orthogonalise after each iteration, 
leaving us with the method of Benettin et al. 31 1. Similar 
observations hold true for the higher Lyapunov spectrum. 

At this point an objection may be raised. Since 



(ao |7773 



t-rt 



.TtT ■ 



•7- 2 Ti|ao), 



(54) 



where the dagger denotes transposition, the correlation 
function (|5^) corresponds to a realisation of the random- 
ness that is always symmetric around the midpoint of 
r\ and r 2 . From the above physical picture leading to 
Eq. it seems that what we really ought to compute 



(b„|b„) 



(55) 



we find that 



where the (transposed) TMs used to obtain the states on 
the left implement a different realisation of the random- 
ness than that used to obtain the states on the right. 

Numerically we are now facing the problem of com- 
puting the average of huge numbers that are no longer 
necessarily positive. As discussed in connection with 
Eq. (43) we do not expect the correlation function to 
be self-averaging, and because of possible negative val- 
ues of Eq. (55) the subterfuge of averaging its logarithm 
will not help us out. Trial runs seem to indicate that 
for sufficiently small values of q and R (such as q = 3 
and R = 2) the matrix elements appearing in Eq. j55| ) 
computed for the usual samples of lcnght 200 may have 
either sign, but that their quotient is invariably posi- 
tive. The corresponding result for xt is roughly equal 
to that obtained from the cumulant expansion. Unfortu- 
nately, for larger values of q or R rare events of negative 
quotients begin to occur, and any attempt of averaging 
Eq. (55) without resorting to logaritms is hampered by 
such large fluctuations as to render the results insuffi- 
ciently accurate at the very best. Computations along 
these lines, though physically appealing, must therefore 
be abandoned on numerical grounds. 
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Yet another possibility of determining at least an ap- 
proximate value of Xt is through the conformal sum rule 
f33| that for an n-fold replicated system reads 



12 l + E-difaJe- 2 " 



(56) 



where the sum runs over all operators in the theory, in- 
cluding the descendants of the Verma module with their 
appropriate multiplicities, and di(n) are the multiplici- 
ties pertaining to the permutational symmetry of the n 
replicas and the q Potts states. For the magnetic opera- 
tor di(n) = n(q — 1), since there are (q — 1) independent 
order parameters, and in the case of the energy operator 
di(n) = n. In the pure system this yields quite accu- 
rate estimates for xt if the exact values of c(l) and xh 
are inserted along with the first descendant of the latter. 
Differentiating and going to the replica limit we find that 
for a random system 
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xi(q- 1) 



-2wxi 



X 2 



(9-1) 



2 c -27rx 2 



(57) 



so that for values of X\, Xi and xt near those of the 
Ising model the term with xt enters only as a small cor- 
rection. Consequently, at the very best only X2 can be 
determined with some confidence from our previous re- 
sults for c'(0) and x\. Its value appears to be consistent 
with that obtained from the cumulant expansion. 

Finally we should like to mention that preliminary 
studies of exact partition function zeros for small L x L 
lattices with quenched bond randomness hints at an in- 
teresting new method of estimating xt- Although the 
different realisations of the quenched bond randomness 
in general lead to a considerable scatter in the positions 
of such zeros, it turns out that the zeros that are closest to 
the real axis only exhibit a very weak dependence on the 
realisation. But these zeros are precisely those that fix 
xt through their finite-size scaling. Results along these 
lines, both for zeros of the Lee- Yang and of the Fisher 
type, will be published elsewhere [p4[. 



F. Phenomenological renormalisation 

In view of the difficulties encountered in our attempts 
to extract xt directly at the critical point we turn our 



attention to the method of phenomenological renormali- 
sation pl[| , which is closer in spirit to the ideas that lead 
to the bound @. 

The magnetic correlation length can be found from the 
TM spectra through 



ti^T)- 1 = In (jjp) = L(fS 2 ~ / n ), (58) 

and we note that this quantity would be self-averaging in 
the random case. Motivated by the form £ ~ (T — T c )~ v 
of the divergence of the correlation length in the infinite 
system we make the finite-size scaling ansatz 



Z{L,t) = L<t>{{T-T tl )L 1 ' v ). 



(59) 



For pure systems, then, one traditionally scans through 
the vicinity of T c to find an effective T C (L) as the solution 
of 



Z(L,T C (L)) _ g(L-l,T c (£)) 
L L-l 

and computes an approximant v(V) from 



1 



1 



In (n{L,T)/»(L-l,T)) 



ln(L/(L-l)) 



(60) 



(61) 



T=T C (L) 



where the derivatives 
PV ' ' dT 



^<j>'({T-T c )L l l v ) (62) 



As L 



oo we 



are found by numerical differentiation, 
have T C (L) — > T c and v{L) — > v. 

In the random case the extracted values of £(L, T) are 
hampered by substatial error bars, and the method just 
outlined becomes by far too inefficient. Fortunately the 
very costly idea of scanning for T C (L) can be discarded, 
since the exact T c of the RBPM is known from Eq. (0). 
Consequently the derivatives (|6^) and the approximants 
( |6l| ) are evaluated at the exact T c , whence the only re- 
maining source of errors is that of statistical fluctuations 
over the different realisations of the randomness. 

Naively one would now find the derivative ( |62| ) by sub- 
traction of the free energies evaluated at T = T c (l ± e), 
where e <C 1. Although this method yields reasonable 
results for e ~ 10~ 2 it is way too inaccurate, since it in- 
volves the subtraction of almost identical quantities (with 
error bars) . A superior strategy is to divide the strip into 
patches of length 200, calculate exact 2 values of /J,(L) 



2 Since we are now faced with differentiating a quantity that 
is known with full machine precision (10 -16 ) we can con- 
centrate on minimising the rounding and truncation errors. 
This is accomplished by choosing e = 10~ . See Sect. 5.7 of: 
W. H. Press et at, Numerical Recipes in C, second edition 
(Cambridge University Press, Cambridge, 1992). 
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for each of those, and finally average over the totality 
of such patches. In this way one exploits the fact that 
£{L,T C (1 — e)) and £(L,T c (l + e)) are strongly correlated 
when the realisation of the randomness is kept fixed. In 
practice we found that this trick lead to a reduction of 
the error bars with a factor ~ 120. 

Sample results obtained by using these prescriptions 
are shown in Table VII . It is seen that although the con- 
vergence is still too slow for reliable extrapolations to the 
limit of an infinite system to be made, the conflict with 
the bound ( fl9| ) appears to be resolved. 

We have found that the convergence of the estimators 
v{L) can be significantly sped up if one performs the 
numerical differentiation ( [62] ) by going perpendicularly 
off the self-duality criterion in (K\,Kz) space instead 
of maintaining the condition K 2 = RK\. Indeed, one 
may imagine that there is another exponent associated 
with motion along the critical surface, and maintaining 
K2 = RK\ one would then measure an admixture of 
this spurious exponent, in particular for large R. A sim- 
ple calculation using Eq. (||) shows that one should then 



evaluate £(L, K[, K' 2 ) at 



K[=K 1 



Re Kl , K 



qe 



A'., 



K! 2 = K 2 {l±e). 



(63) 



The sample results shown in Table VIII exhi bit a con- 
spicuous improvement over those of Table VII . Not only 
is the convergence faster, but it is even seen that the es- 
timators v(V) form a monotonically increasing sequence 
for low values of R and a monotonically decreasing one 
for high R. The extrapolated v is pinched between those 
two sequences and consequently quite accurately deter- 
mined. 

Plots of the estimators v(V) for 3 < L < 7 and several 
values of R are shown in Figs. || and |?] for q = 8 and 
q = 64 respectively. It is seen that the curves for v(V) 
and v(L — 1) intersect at a unique value of R, seemingly 
converging quite fast to a definite value R* as L increases. 
We interpret R* as the randomness strength at the crit- 
ical FP and the corresponding value of v as the correct 
critical exponent. It is tempting to conjecture that the 
curves v{L) approach v on the entire interval R €]l,oo[ 
as L — > 00. From the graphs it seems that the conver- 
gence is faster for large q. 

The values of v and R* corresponding to this scenario 



are shown in Table IX . In accordance with the phase dia- 
gram (Fig. |l|) R* is a slowly, supposedly logarithmically 3 
, increasing function of q. For q = 2 the deviation from 



v = 1 can be ascribed to logarithmic corrections Hj , and 
for q = 3 our result for v is in agreement with the (q — 2)- 
expansion pj] though the possibility of replica symmetry 
breaking cannot be ruled out . Also for q > 4 are our 
values for v numerically consistent with unity, indicating 
that, unlike what is the case for the magnetic expononent, 
the thermal one displays only a weak (/-dependence. 

From Figs. |6| and |^ a remarkable feature about the pure 
system (R = 1) is apparent. For q = 8 the estimators 
v{L) seem to converge to v = | whilst for q = 64 the 
extrapolated value is v ~ 0. The former value is hardly 
surprising since, as we also remarked above, a first-order 
transition is expected to exhibit scaling with trivial ef- 
fective exponents (in this case: xt = 0). On the other 
hand, v ~ for q = 64 has to do with the length scales 
of the system. Namely, from the asymptotic behaviour 
of the correlation length of the pure system Mw 



Inq 



as q 



(64) 



we infer that £ ~ 1 -C L at the transition point of the 
q = 64 system. But away from the transition point we 
also expect £ = 0(1), since the lattice spacing (unity) 
is the least length scale in the system. After all there 
is a ferromagnetic interaction between nearest neighbour 
spins. We thus conclude that £ is roughly temperature 
independent in this case. In order for this to be con- 
sistent with the asymptotic behaviour of the finite-size 
scaling ansatz ([sol) 



4>{x) 



~ v for x < 1 



(65) 



we must then have v ~ 0. This is to be contrasted to the 
case of q = 8 where ^ > L so that we clearly "see" the 
phase transitions in our strips of width L. 

Very recently the bound (^) was challenged by 
Pazmandi et al. p4| . These authors claimed that the 
standard method of averaging over the disorder in finite- 
size (FS) systems introduces a new diverging length scale 
into the problem, whence the resulting j/ps may be unre- 
lated to the true exponent v governing the divergence of 
the correlation length in the infinite system. In particular 
v can be less than 1 , and if this is the case the standard 
method is liable to yield exactly z^pg = |. Ref. |24| then 
went on to suggest a noise reduction procedure that pro- 
fessedly would allow one to access the true v. For each 
realisation of the binary randomness (^) used in the dis- 
order average there will be a fraction p of weak bonds 
(K\ ) . The noise due to the fluctuations of p around its 
average value V = \ can then be reduced by adjusting the 



3 This supposition constitutes the simplest possibility al- 
lowed by the phase diagram of Fig. [I] in which R* — ► 00 as 
lng — > 00. 
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couplings (Ki,K 2 ) for that particular realisation to the 
values they would assume at the critical point of an infi- 
nite system with a (mean) fraction of weak bonds equal 
to p. 

To implement this we are faced with the task of find- 
ing the two-dimensional critical surface in the space 
(Ki , K 2 , p) using our knowledge of its one-dimensional 
intersection with the plane p = 5, viz. Eq. (^|). Let the 
fraction of weak bonds in a particular realisation be 



5(1 + Cp). 



(66) 



where e p <C 1. The symmetry p <-> 1 — p ensures that, to 
first order in e p , we must still go perpendicularly off the 



self-dual curve in the (Ki,K 2 ) subspace, as in Eq. (63). 
Since an increase in the number of weak bonds must be 
offset by an increase of the K' s in order to keep the cou- 
pling to the energy density constant, the correct prescrip- 
tion is 



K x -> K x \\ + e K 
K 2 -» K 2 (l + e K ) 



Re R i 



ge 



K 2 



(67) 



for some ex > 0. Demanding that the combined change 
in p and in (Ki,K 2 ) must leave the coupling to the en- 
ergy density invariant furnishes a relation between e p and 



K 2 - K x 



P ts «ef , K 



l) 2 + K 2 



(68) 



and the derivatives ( p2| ) are now evaluated at these val- 
ues of the parameters by going perpendicularly off the 
critical surface. To first order, of course, Eq. ( |63|) still 
gives the correct way of doing so. 



Our confidence in the results of Table |IX| is increased by 
observing that the implementation of this novel averaging 
procedure does not alter our results. Indeed, trial runs for 
q = 8, where the discrepancy between the xt extracted 
from the Lyapunov spectrum and phenomenological RG 
respectively is large, render the values of the estimators 
/x(L) unchanged within (small) error bars. It is thus con- 
cluded that even though our results for v are conspic- 
uously close to satisfying the bound (E|) with equality, 
this is not due to an artifact in the averaging procedure. 



G. The higher Lyapunov spectrum 



the first five gaps of the i^-even sector are related to the 
energy density e, its first descendants L_ie and L_ie, 
the stress tensor T and its conjugate T. To wit, the scal- 
ing dimensions of these operators can be found from the 
gaps through Eq. (Q) , and we have verified this using our 
connectivity basis TMs. 

In view of the bound (|4^) it is problematic to associate 
the first gap with the energy density in the random case, 
but it is nevertheless a beguiling question whether such 
concepts as descendants and the stress tensor are pre- 
served by the randomness. To investigate this issue we 
have computed the first few Lyapunov exponents of T 11 
for 1 < L < 8, averaging over 100 runs as usual. The 
scaling dimensions corresponding to the first five gaps 
for q = 3,4,5 and R — 2 are shown in Table [j§ Self- 
averaging was ensured by utilising the cumulant expan- 
sion, and parabolic least-squares fits included the first 
three cumulants. 

It is quite remarkable that even if the scaling dimen- 
sion corresponding to the first gap may not be equal to 
xt our data give strong reasons to believe that it has a 
descendant, and that this descendant has the expected 
degeneracy of two. And even though the scaling dimen- 
sions in general depend on q those corresponding to the 
fourth and the fifth gaps are constant within error bars 
and very close to 2, as is expected for the stress ten- 
sor of a conformally invariant system p3| . Preliminary 
data for even higher Lyapunov exponents seem to hint at 
descendants at level two, but since we have found that 
in the pure system higher and higher eigenvalues require 
larger and larger system sizes before the asymptotic scal- 
ing form (^) is valid, massive computations are needed to 
establish reliable results for all but the first few scaling 
dimensions. 

Another interesting feature of our data for the higher 
Lyapunov exponents is that the Harris criterion seems to 
be valid in a very complete sense. Namely, trial runs for 
q = 1 seem to indicate that although individual cumu- 
lants exhibit a pronounced dependence of R, their sum is 
virtually independent of the strength of the randomness 
in the whole range R G [1,2]. It thus appears that all 
exponents Xi that we can extract numerically from the 
Lyapunov spectrum, using Eq. (Q) and the cumulant ex- 
pansion, obey the Harris criterion. Since the connection 
between these exponents and the scaling dimensions of 
the underlying CFT is not completely known (witness 
xt) this may well turn out to be a non-trivial observa- 
tion. 



Although the second Lyapunov exponent of T 11 fails 
to yield the thermal scaling dimension xt in the standard 
way it is hard to believe that the Lyapunov spectrum is 
not in some way related to the operator content of the 
CFT underlying the RBPM. In the case of the pure three- 
state Potts model, for example, it is well known Ga] that 



V. DISCUSSION AND OUTLOOK 

In a recent paper by Picco E7J] it has been suggested 
that for q — TP the effective central charge at the ran- 
dom FP is d = §, and that this class of models thus 
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behaves as p decoupled Ising models. Even without re- 
ferring to our values of the magnetic exponent we should 
like to point out that all the data show is that d oc In q 
with a constant of proportionality that is very close to 
2)^-2 — 0.721. But this constant is also very close to that 

of the percolation point, viz. ~ 0.689. Indeed, these 
two numbers differ by less than 5%, and since our error 
bars and those of Picco are in the 2% and the 4% range 
respectively, there is no irrefutable way of numerically 
distinguishing between percolative, Ising-like or indeed 
some other, presently unknown, behaviour of the central 
charge. A similar observation is valid for 2 < q < 4 where 
our numerical data as displayed in Fig. [| are compatible, 
within error bars, with both the values at the pure and 
the random FP (but not, in this case, with those at the 
percolation point). 

On the other hand, our results for the magnetic expo- 
nent should leave no doubt that the correct CFT describ- 
ing the RBPM cannot be that of a number of decoupled 
Ising models. In particular, the non-Ising value at q = 8 
is in sharp disagreement with the Monte Carlo results of 
Rcf . . One possible explanation of this discrepancy is 
that these authors define a non-standard order parameter 
through 

p = L'- d (m^(N 1 ,N 2 ,...,N q )), (69) 

where Ni is the number of spins in state i, which is re- 
lated to our local order parameter defined in Eq. ( p6|) by 
Ni = ^2 r (Mi(r) + q^ 1 ). The site label r runs over a hy- 
percube of side L' with 24 < L' < 84. Note that p may 
also be written as 

P = L^ &t (±(N} ) y . (70) 

Expressed in terms of the local order parameter, (Nf) 
gives a sum of terms each of the form 

Y.^Mri^M^p . ..Mi(r„) fc -), (71) 

where fci + k% + • • • = k. As k — > oo at fixed L', it is 
clear that at least some of the kj must grow large. In 
the pure system, this should not matter, since each term 
will scale in the same manner. But when multiscaling is 
present, the scaling behaviour of the kj power of the lo- 
cal order parameter may be different. Indeed, in the limit 
of k — > oo one would expect p to scale with dimension 
limbec Xk/k, which is less than xi by convexity. 

Another criticism of Ref. |[l| is that the realisations of 
the binary randomness considered were confined to those 
for which the number of strong and weak bonds in each 
of the two lattice directions were equal. Though this 
restriction is clearly innocuous in the limit L' —> oo this 
may not be so as far as the finite-size scaling is concerned. 



From trial runs where similar restrictions were imposed 
to the bond distributions of the TMs we found that seem- 
ingly harmless noise reductions schemes can influence the 
output substantially. 

Finally, the mapping to the RFIM p3| illustrated that 
for large q typical configurations consist of large clusters 
of spins in the same Potts state, separated by domain 
walls. Whilst our very long strips are guaranteed to ac- 
comodate large regions in which all q values of the order 
parameters are realised, it is not clear that this should 
be the case in the much smaller square geometries of 
Ref. |ll[] . Indeed it seems likely that one would find Ising 
exponents if the geometry under consideration typically 
can accomodate at most two different large clusters. 

We now turn our attention to the thermal exponent. If 
the phenomenological RG scheme is to be trusted the val- 
ues of xt only exhibit a weak dependence on q, although 
the (q— 2)-expansion gives us reason to believe that there 
is some variation ||l5| . It is interesting that xt stays so 
close to unity even at very high q, but presently we do 
not have any arguments to explain this finding. Unfor- 
tunately the method employed was unable to resolve the 
slight deviations from unity, and it is indeed a challenge 
to future research to find more accurate ways of assess- 
ing xt for disordered systems. Our results on the higher 
Lyapunov spectrum are nothing if not intriguing, and we 
believe that a great effort must be made to understand 
why the first gap in the spectrum fails to be related to xt 
in the standard way, even though the higher gaps show 
clear indications of a conformal field theory underlying 
the RBPM. 

A very interesting issue to be addressed by future 
research is that of the dynamical universality class of 
the RBPM. In particular it would be interesting to see 
whether the dynamic critical exponent z does or does 
not agree with the Ising value of z m 2, or whether, in 
analogy with the RFIM, there is logarithmically slow dy- 
namics due to the pinning of domain walls by impurities. 

Other types of randomness are also of interest to the 
question whether a first-order phase transition is softened 
due to impurities. In this paper we have studied the effect 
of quenched bond randomness in a flat, regular lattice. A 
somewhat different scenario is obtained by investigating 
the pure g-state Potts model on lattices with quenched 
connectivity disorder. In Ref. ]5^ | MC simulations of the 
q = 8 model on two-dimensional Poissonian random lat- 
tices (Voronoi tessellations) with toroidal topology un- 
ambiguously showed that the first-order nature of the 
transition was not modified. 

An argument why this must be so, at least for large q, 
can readily be given. For simplicity consider the model 
on the dual Delaunay random lattice, which per con- 
struction is a triangulation of space p6[ |. As on the 
regular lattice, at large q there are only two important 
states in the equivalent random cluster model: the empty 
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lattice, which contributes a term ^vertices to the par- 
tition function, and the full lattice, contributing a fac- 



tor u Ncdgos , where u — e K — 1. Since for any triangu- 
lation 2iV dgos = 3-/V vcr ti CC s, the transition occurs when 
u ~ g 2 / 3 for any geometry. If we now consider that 
part of the lattice within a large hypercube of side L, 
the fluctuations in the difference of the energies of these 
two states inside this region will come solely from the 
edges which penetrate the boundary. On average, the 
difference in the energies of these two states will still be 
zero, but there will be fluctuations of the order of the 
square root of the number of bonds which penetrate the 
boundary, which will therefore be C(L' d_1 )/ 2 ). These 
are very much smaller than the analogous fluctuations 
which are present when random bonds are added: these 
are 0(L d ^, which leads to the Imry-Ma type of argu- 
ment HHH. For d = 2 Voronoi tessellations the fluc- 
tuations are thus always smaller than the domain wall 
energy C(L d_1 ), and we conclude that such randomness 
is strongly irrelevant (at least for large q), in agreement 
with the results of Ref . |55| . 

Yet another kind of randomness is obtained by study- 
ing the Potts model on quenched random gravity graphs, 
for which MC simulations for q = 10 have provided strong 
evidence for a softening scenario similar to ours [ |57| . 
However, in this case the curvature is random and hence 
when the lattice is embedded in the plane, it is fractal. 
Although the argument about compensation of the bulk 
energies when u ~ q 2 ^ 3 works for any triangulation, the 
number of boundary edges may well scale in a different 
manner for these lattices. Whilst it would be interesting 
to study this in detail, it is clear that this problem is 
quite different from ours, and neither our arguments nor 
those of Refs. |||| can be directly applied. 
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FIG. 1. Schematic phase diagram in the critical surface for d > 2. q increases to the left and w is the disorder strength, 
with P1P2 being the percolation limit. RG flows are indicated. The latent heat is non-vanishing within the shaded region, and 
elsewhere the transition is continuous, controlled by the line of fixed points Piqi. As d — + 2 the shaded region collapses to a 
line q^O of first-order transitions in the pure system. 

FIG. 2. Plots of X(L), normalised to A(l) = 1, showing that bond randomness renders the phase transition second order. 
The random systems have R = K2/K1 — 2. 

FIG. 3. The effective central charge c' as a function of q for 2 < q < 4. The perturbative results by Ludwig and Cardy [14] 
do not differ appreciately within the range of q- values where the expansion is supposed to be valid. Accordingly the numerical 
data are unable to distinguish between pure and non-trivial random behaviour. They are also quite close to, but distinguishable 
from, the percolation point values. 

FIG. 4. Effective central charge as a function of log 10 q for 1.5 < q < 64. For large q the data for R — 10 are supposed to 
represent the true behaviour at the random fixed point, as argued in the text. 

FIG. 5. Magnetic exponent xi = f3/v for R = 2 as obtained from the cumulant expansion. x\ is an increasing function of q, 
continuously connecting onto the perturbative results near q = 2. For q > 8 the cumulant expansion begins to break down. 

FIG. 6. Estimants v(V) for the thermal exponent at q — 8 as obtained from phenomenological renormalisation applied to 
strips of width L and L — 1. In the pure system (R — > 1, see rightmost inset) the estimants converges to 1/(00) = |. Curves for 
neighbouring system sizes intersect at values of v and R that converge to those at the random fixed point as L — * 00. In this 
case v — 1.01 ± 0.02 and R* » 9 (see left inset). Error bars are no larger than the size of the symbols. 

FIG. 7. Phenomenological renormalisation at q — 64. The curves intersect at larger angles than before, allowing for a rather 
accurate determination v = 1.02 ±0.03 in spite of the large fluctuations. Error bars are comparable to the size of the symbols. 
From the rightmost inset it is seen that v — > in the pure systems, as explained in the text. The left inset is a magnification 
of the region around R* ^ 10. 
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2,496,144 


4,213,597 



TABLE I. The number of connectivity states for a Potts model transfer matrix of width L with (di) and without (cl) an 
external magnetic field. Also shown is the size of the magnetic sector when using a ghost site (cLl — cl) and a seam (Lcl)- For 
large strip widths the seam is advantageous. The number cl of well-nested L-point connectivities should be compared to the 
total number of L-point connectivities bL which increases faster than exponentially as a function of L. 
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(8) 


3.48241 


(9) 


4.67423 


(11) 


4 


1 


.89550 


(6) 


2. 


17203 


(6) 


2 


,37431 


(6) 


2.88328 


(7) 


3.42387 


(8) 


4.59557 


(10) 


5 


1 


.88895 


(5) 


2. 


16182 


(6) 


2 


,36126 


(6) 


2.86392 


(6) 


3.39934 


(7) 


4.56442 


(9) 


6 


1 


.88568 


(5) 


2. 


.15649 


(5) 


2 


,35449 


(5) 


2.85395 


(6) 


3.38683 


(6) 


4.54893 


(8) 


7 


1 


.88377 


(4) 


2. 


.15328 


(5) 


2 


,35040 


(5) 


2.84798 


(5) 


3.37948 


(6) 


4.53974 


(7) 


8 


1 


.88250 


(4) 


2. 


.15113 


(4) 


2 


,34782 


(4) 


2.84424 


(5) 


3.37479 


(5) 


4.53394 


(7) 


9 


1 


.88164 


(4) 


2. 


14993 


(24) 


2 


,34624 


(25) 


2.84172 


(11) 


3.37186 


(31) 


4.53017 


(39) 


10 


1 


.88098 


(4) 


2. 


.14858 


(23) 


2 


,34504 


(23) 


2.84011 


(26) 


3.36918 


(29) 


4.52653 


(25) 


11 


1 


.88048 


(4) 


2. 


14804 


(22) 


2 


,34386 


(22) 


2.83851 


(24) 


3.36768 


(29) 


4.52465 


(34) 


12 


1 


.88017 


(3) 


2. 


14744 


(20) 


2 


,34314 


(21) 


2.83765 


(24) 


3.36639 


(26) 


4.52316 


(34) 


13 


1 


.87991 


(3) 



























TABLE II. Critical free energies per site, — /"(L), for R — 2 and various values of q. The figures in parentheses indicate 
the error bar on the last quoted digits. 



L 


g = 2 


9 = 3 


q = 4 


9 = 8 


q = 8 


q = 16 


q = 16 


q = 64 


q = 64 




R = 2 


R = 2 


R = 2 


7? = 2 


R=W 


R = 2 


R= 10 


R=2 


R= 10 


1 


0.563 (1) 


0.927 (1) 


1.184 (1) 


1.787 (1) 


1.731 (3) 


2.330 (1) 


2.322 (4) 


3.120 (1) 


3.476 (5) 


2 


0.508 (2) 


0.825 (3) 


1.042 (2) 


1.515 (3) 


1.586 (10) 


1.864 (3) 


2.101 (10) 


2.194 (4) 


3.150 (13) 


3 


0.499 (3) 


0.800 (6) 


1.003 (6) 


1.441 (7) 


1.521 (23) 


1.752 (8) 


2.052 (25) 


2.065 (10) 


3.034 (30) 


4 


0.500 (6) 


0.813 (14) 


0.994 (14) 


1.424 (15) 


1.548 (52) 


1.750 (17) 


2.089 (57) 


2.157 (22) 


3.079 (68) 


5 


0.505 (11) 


0.842 (30) 


1.005 (31) 


1.426 (33) 


1.622 (113) 


1.785 (38) 


2.203 (125) 


2.305 (47) 


3.209 (148) 


6 


0.500 (20) 


0.818 (62) 


0.963 (63) 


1.360 (67) 


1.587 (228) 


1.794 (78) 


2.196 (251) 


2.384 (93) 


3.213 (300) 



TABLE III. Effective central charge d extracted from parabolic fits with Lo < L < L max , as described in the text. Error 
bars on the last quoted digit are shown in parentheses. The choice Lo = 3 appears to be optimal, provided that the strength 
of the randomness R is large enough (see text), and the corresponding values of c , shown in bold face, should be regarded as 
our results. 
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P 9 = 2 


g = 4 


q = 8 


q = 16 


q = 32 


q = 64 


0.05 0.522 (25) 


1.030 (26) 


1.477 (27) 


1.780 (29) 


1.920 (30) 


1.974 (31) 


0.10 0.519 (35) 


1.032 (36) 


1.510 (38) 


1.915 (40) 


2.251 (42) 


2.552 (43) 


0.15 




1.539 (46) 


1.996 (48) 


2.416 (50) 




TABLE IV. Effective central charge c obtained using a trinary distribution of the random bonds. There is 


a fraction p of 


very weak and very strong bonds respectively, the remaining fraction 1 


— 2p being critical. 






L g = 2 


q = 3 


q = 4 


q = 8 


q = 16 


~ P. A 

q = 04 


R = 2 


R = 2 


R = 2 


R=10 


R=W 


p in 


2 0.5662 (6) 


0.9300 (7) 


1.1903 (7) 


1.723 (3) 


2.313 (3) 


9 A RCl I A \ 
0.409 (4J 


3 0.535 (1) 


0.878 (1) 


1.117 (1) 


1.656 (5) 


2.207 (6) 


O Q 1 1 ( r 7\ 

o.oll ((J 


4 0.521 (2) 


0.848 (3) 


1.075 (3) 


1.605 (10) 


2.148 (11) 


o one; ( i o\ 
o.zU6 (16) 


5 0.514 (3) 


0.836 (5) 


1.051 (5) 


1.585 (19) 


2.125 (20) 


3.163 (24) 


6 0.512 (4) 


0.839 (9) 


1.043 (9) 


1.595 (34) 


2.144 (38) 


3.174 (45) 


7 0.510 (6) 


0.840 (18) 


1.022 (18) 








8 0.509 (9) 


0.812 (33) 


1.011 (34) 








9 0.503 (14) 












10 0.500 (23) 













TABLE V. Effective central charge c' extracted from linear fits of /^(L) — /"(co) versus l/L 2 , with Lo < L < I/ max . In all 
cases the approach towards the asymptotic values of Table III is from above. Error bars on the last quoted digit are shown in 
parentheses. 





L 


1. cumulant 


2. cumulant 


3. cumulant 


4. cumulant 


5. cumulant 


6. cumulant 


9 = 3 


1 


-1.039786 (242) 


0.060716 


-0.000791 


0.000830 


-0.004583 


-0.000778 




2 


-0.253209 (146) 


0.012391 


-0.000347 


0.000279 


0.000059 


0.000386 




3 


-0.106163 (113) 


0.004963 


-0.000246 


0.000063 


-0.000102 


0.000046 




4 


-0.057901 (95) 


0.002784 


-0.000143 


0.000006 


0.000034 


0.000003 




5 


-0.036521 (84) 


0.001810 


-0.000105 


0.000001 


-0.000002 


-0.000003 




6 


-0.025172 (76) 


0.001289 


-0.000075 


0.000008 


-0.000001 


-0.000002 






-0.018426 (69) 


0.000968 


-0.000069 


0.000002 


0.000002 


-0.000001 


q = 8 


1 


-1.380171 (289) 


0.104382 


0.004069 


0.014001 


0.019452 


-0.013889 




2 


-0.326484 (177) 


0.028366 


-0.001683 


-0.000432 


0.000157 


-0.003145 




3 


-0.132560 (138) 


0.014908 


-0.001822 


0.000356 


0.000221 


-0.000083 




4 


-0.071296 (115) 


0.010129 


-0.001610 


0.000319 


-0.000959 


0.002323 




5 


-0.044886 (102) 


0.007880 


-0.001619 


0.000252 


-0.000082 


-0.000456 




6 


-0.031045 (92) 


0.006450 


-0.001538 


0.000607 


-0.000184 


-0.001096 




7 


-0.022851 (84) 


0.005401 


-0.001505 


0.000237 


0.000234 


-0.000280 



TABLE VI. The first six cumulants of — A/(L) for 1 < L < 7 and R — 2. The error bar on the first cumulant (shown in 
parentheses) is related to the second cumulant; error bars on the higher cumulants are not shown. 
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L /z(L) v(L) 

1 1.087 (1) 

2 4.229 (3) 1.041 (1) 

3 10.426 (8) 0.816 (2) 

4 19.682 (18) 0.827 (3) 

5 31.867 (33) 0.863 (5) 

6 46.994 (53) 0.885 (7) 

7 65.020 (79) 0.904 (9) 



TABLE VII. Phenomenological renormalisation for the thermal scaling dimension xt — 2 — \/v at q = 8 and R — 10. For 
each strip width L the 100 independent strips of length m = 10 5 are divided into patches of length 200. Within each patch the 
exact jj,(L,T c ) is computed, based on evaluations of £(L, K' 1: K'2) at K[ = Ki(l ± €k) and K' 2 = RK[, where K\ is found from 
Eq. (3). Final results and error bars are obtained as the mean value and the standard deviation over the totality of patches. 





L 


lx(L) 


u(L) 


R = 6 


2 


1.898 (1) 






3 


4.456 (2) 


0.905 (1) 




4 


8.172 (6) 


0.902 (2) 




5 


12.974 (11) 


0.933 (4) 




6 


18.883 (18) 


0.945 (6) 




7 


25.891 (27) 


0.955 (8) 


R = 10 


2 


1.832 (1) 






3 


3.917 (2) 


1.144 (2) 




4 


6.815 (5) 


1.081 (4) 




5 


10.486 (9) 


1.074 (6) 




6 


14.948 (15) 


1.059 (8) 




7 


20.198 (22) 


1.050 (11) 



TABLE VIII. Phenomenological renormalisation going perpendicularly off the critical surface for q — 8 and R = 6 and 10 
respectively. The data collection was done as before. 



1 


V 


R* 


2 


1.12 (3) 


7(1) 


3 


1.04 (4) 


8(1) 


8 


1.01 (2) 


9(1) 


64 


1.02 (3) 


10 (1) 



TABLE IX. Values of the critical exponent v and the randomness strength R* at the random fixed point as obtained from 
phenomenological renormalisation. 



q 1. gap 2. gap 3. gap 4. gap 5. gap 



3 


0.899 (4) 


1.877 (13) 


1.885 (12) 


2.045 (24) 


2.050 (23) 


4 


0.817 (5) 


1.811 (9) 


1.818 (8) 


2.043 (23) 


2.049 (24) 


5 


0.754 (6) 


1.771 (6) 


1.779 (6) 


2.058 (24) 


2.065 (25) 



TABLE X. Scaling dimensions corresponding to the first five gaps in the Lyapunov spectrum of T for R = 2. The parabolic 
least-squares fits included the first three cumulants of the probability distribution, and error bars were extracted based on the 
fits with Lo = 4, 5 and 6. 
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